MOAB: Mesh Oriented datABase  (version 5.4.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,
00342                                         MeshDomain* domain,
00343                                         OFEvaluator& obj_eval,
00344                                         const Settings* settings,
00345                                         MsqError& err )
00346 {
00347     const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
00348     PatchData global_patch;
00349     if( settings ) global_patch.attach_settings( settings );
00350 
00351     // if we need to fill out the global patch data object.
00352     if( ( totalFlag & ( GRAD_FLAGS | OF_FLAGS | VERTEX_MOVEMENT_RELATIVE | UNTANGLED_MESH ) ) || timeStepFileType )
00353     {
00354         global_patch.set_mesh( mesh );
00355         global_patch.set_domain( domain );
00356         global_patch.fill_global_patch( err );MSQ_ERRRTN( err );
00357     }
00358 
00359     // now call the other reset
00360     reset_inner( global_patch, obj_eval, err );MSQ_ERRRTN( err );
00361 }
00362 
00363 /*!Reset function using using a PatchData object.  This function is
00364   called for the inner-stopping criterion directly from the
00365   loop over mesh function in VertexMover.  For outer criterion,
00366   it is called from the reset function which takes a MeshSet object.
00367   This function prepares the object to be used by setting the initial
00368   values of some of the data members.  As examples, if needed, it resets
00369   the cpu timer to zero, the iteration counter to zero, and the
00370   initial and previous objective function values to the current
00371   objective function value for this patch.
00372   The return value for this function is similar to that of terminate().
00373   The function returns false if the checked criteria have not been
00374   satisfied, and true if they have been.  reset() only checks the
00375   GRADIENT_INF_NORM_ABSOLUTE, GRADIENT_L2_NORM_ABSOLUTE, and the
00376   QUALITY_IMPROVEMENT_ABSOLUTE criteria.  Checking these criteria
00377   allows the QualityImprover to skip the entire optimization if
00378   the initial mesh satisfies the appropriate conditions.
00379  */
00380 void TerminationCriterion::reset_inner( PatchData& pd, OFEvaluator& obj_eval, MsqError& err )
00381 {
00382     const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
00383 
00384     // clear flag for BOUNDED_VERTEX_MOVEMENT
00385     vertexMovementExceedsBound = 0;
00386 
00387     // Use -1 to denote that this isn't initialized yet.
00388     // As all valid values must be >= 0.0, a negative
00389     // value indicates that it is uninitialized and is
00390     // always less than any valid value.
00391     maxSquaredMovement = -1;
00392 
00393     // Clear the iteration count.
00394     iterationCounter = 0;
00395 
00396     // reset the inner timer if needed
00397     if( totalFlag & CPU_TIME )
00398     {
00399         mTimer.reset();
00400     }
00401 
00402     // GRADIENT
00403     currentGradInfNorm = initialGradInfNorm = 0.0;
00404     currentGradL2NormSquared = initialGradL2NormSquared = 0.0;
00405     if( totalFlag & GRAD_FLAGS )
00406     {
00407         if( !obj_eval.have_objective_function() )
00408         {
00409             MSQ_SETERR( err )
00410             ( "Error termination criteria set which uses objective "
00411               "functions, but no objective function is available.",
00412               MsqError::INVALID_STATE );
00413             return;
00414         }
00415         int num_vertices = pd.num_free_vertices();
00416         mGrad.resize( num_vertices );
00417 
00418         // get gradient and make sure it is valid
00419         bool b = obj_eval.evaluate( pd, currentOFValue, mGrad, err );MSQ_ERRRTN( err );
00420         if( !b )
00421         {
00422             MSQ_SETERR( err )
00423             ( "Initial patch is invalid for gradient computation.", MsqError::INVALID_STATE );
00424             return;
00425         }
00426 
00427         // get the gradient norms
00428         if( totalFlag & ( GRADIENT_INF_NORM_ABSOLUTE | GRADIENT_INF_NORM_RELATIVE ) )
00429         {
00430             currentGradInfNorm = initialGradInfNorm = Linf( mGrad );
00431             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Initial gradient Inf norm: "
00432                                              << " " << RPM( initialGradInfNorm ) << std::endl;
00433         }
00434 
00435         if( totalFlag & ( GRADIENT_L2_NORM_ABSOLUTE | GRADIENT_L2_NORM_RELATIVE ) )
00436         {
00437             currentGradL2NormSquared = initialGradL2NormSquared = length_squared( mGrad );
00438             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Initial gradient L2 norm: "
00439                                              << " " << RPM( std::sqrt( initialGradL2NormSquared ) ) << std::endl;
00440         }
00441 
00442         // the OFvalue comes for free, so save it
00443         previousOFValue = currentOFValue;
00444         initialOFValue  = currentOFValue;
00445     }
00446     // find the initial objective function value if needed and not already
00447     // computed.  If we needed the gradient, we have the OF value for free.
00448     // Also, if possible, get initial OF value if writing plot file.  Solvers
00449     // often supply the OF value for subsequent iterations so by calculating
00450     // the initial value we can generate OF value plots.
00451     else if( ( totalFlag & OF_FLAGS ) ||
00452              ( plotFile.is_open() && pd.num_free_vertices() && obj_eval.have_objective_function() ) )
00453     {
00454         // ensure the obj_ptr is not null
00455         if( !obj_eval.have_objective_function() )
00456         {
00457             MSQ_SETERR( err )
00458             ( "Error termination criteria set which uses objective "
00459               "functions, but no objective function is available.",
00460               MsqError::INVALID_STATE );
00461             return;
00462         }
00463 
00464         bool b = obj_eval.evaluate( pd, currentOFValue, err );MSQ_ERRRTN( err );
00465         if( !b )
00466         {
00467             MSQ_SETERR( err )
00468             ( "Initial patch is invalid for evaluation.", MsqError::INVALID_STATE );
00469             return;
00470         }
00471         // std::cout<<"\nReseting initial of value = "<<initialOFValue;
00472         previousOFValue = currentOFValue;
00473         initialOFValue  = currentOFValue;
00474     }
00475 
00476     if( totalFlag & ( GRAD_FLAGS | OF_FLAGS ) )
00477         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Initial OF value: "
00478                                          << " " << RPM( initialOFValue ) << std::endl;
00479 
00480     // Store current vertex locations now, because we'll
00481     // need them later to compare the current movement with.
00482     if( totalFlag & VERTEX_MOVEMENT_RELATIVE )
00483     {
00484         if( initialVerticesMemento )
00485         {
00486             pd.recreate_vertices_memento( initialVerticesMemento, err );
00487         }
00488         else
00489         {
00490             initialVerticesMemento = pd.create_vertices_memento( err );
00491         }
00492         MSQ_ERRRTN( err );
00493         maxSquaredInitialMovement = DBL_MAX;
00494     }
00495     else
00496     {
00497         maxSquaredInitialMovement = 0;
00498     }
00499 
00500     if( terminationCriterionFlag & UNTANGLED_MESH )
00501     {
00502         globalInvertedCount = count_inverted( pd, err );
00503         // if (innerOuterType==TYPE_OUTER) MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o
00504         // Num Inverted: " << " " << globalInvertedCount << std::endl;
00505         patchInvertedCount = 0;MSQ_ERRRTN( err );
00506     }
00507 
00508     if( timeStepFileType )
00509     {
00510         // If didn't already calculate gradient abive, calculate it now.
00511         if( !( totalFlag & GRAD_FLAGS ) )
00512         {
00513             mGrad.resize( pd.num_free_vertices() );
00514             obj_eval.evaluate( pd, currentOFValue, mGrad, err );
00515             err.clear();
00516         }
00517         write_timestep( pd, mGrad.empty() ? 0 : arrptr( mGrad ), err );
00518     }
00519 
00520     if( plotFile.is_open() )
00521     {
00522         // two newlines so GNU plot knows that we are starting a new data set
00523         plotFile << std::endl << std::endl;
00524         // write column headings as comment in data file
00525         plotFile << "#Iter\tCPU\tObjFunc\tGradL2\tGradInf\tMovement\tInverted" << std::endl;
00526         // write initial values
00527         plotFile << 0 << '\t' << mTimer.since_birth() << '\t' << initialOFValue << '\t'
00528                  << std::sqrt( currentGradL2NormSquared ) << '\t' << currentGradInfNorm << '\t' << 0.0 << '\t'
00529                  << globalInvertedCount << std::endl;
00530     }
00531 }
00532 
00533 void TerminationCriterion::reset_patch( PatchData& pd, MsqError& err )
00534 {
00535     const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
00536     if( totalFlag & MOVEMENT_FLAGS )
00537     {
00538         if( previousVerticesMemento )
00539             pd.recreate_vertices_memento( previousVerticesMemento, err );
00540         else
00541             previousVerticesMemento = pd.create_vertices_memento( err );MSQ_ERRRTN( err );
00542     }
00543 
00544     if( totalFlag & UNTANGLED_MESH )
00545     {
00546         patchInvertedCount = count_inverted( pd, err );
00547         // MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Patch Inverted: " << " " <<
00548         // patchInvertedCount << std::endl;MSQ_ERRRTN( err );
00549     }
00550 }
00551 
00552 void TerminationCriterion::accumulate_inner( PatchData& pd, OFEvaluator& of_eval, MsqError& err )
00553 {
00554     double of_value = 0;
00555 
00556     if( terminationCriterionFlag & GRAD_FLAGS )
00557     {
00558         mGrad.resize( pd.num_free_vertices() );
00559         bool b = of_eval.evaluate( pd, of_value, mGrad, err );MSQ_ERRRTN( err );
00560         if( !b )
00561         {
00562             MSQ_SETERR( err )
00563             ( "Initial patch is invalid for gradient compuation.", MsqError::INVALID_MESH );
00564             return;
00565         }
00566     }
00567     else if( terminationCriterionFlag & OF_FLAGS )
00568     {
00569         bool b = of_eval.evaluate( pd, of_value, err );MSQ_ERRRTN( err );
00570         if( !b )
00571         {
00572             MSQ_SETERR( err )
00573             ( "Invalid patch passed to TerminationCriterion.", MsqError::INVALID_MESH );
00574             return;
00575         }
00576     }
00577 
00578     accumulate_inner( pd, of_value, mGrad.empty() ? 0 : arrptr( mGrad ), err );MSQ_CHKERR( err );
00579 }
00580 
00581 void TerminationCriterion::accumulate_inner( PatchData& pd, double of_value, Vector3D* grad_array, MsqError& err )
00582 {
00583     // if terminating on the norm of the gradient
00584     // currentGradL2NormSquared = HUGE_VAL;
00585     if( terminationCriterionFlag & ( GRADIENT_L2_NORM_ABSOLUTE | GRADIENT_L2_NORM_RELATIVE ) )
00586     {
00587         currentGradL2NormSquared = length_squared( grad_array, pd.num_free_vertices() );  // get the L2 norm
00588         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- gradient L2 norm: "
00589                                          << " " << RPM( std::sqrt( currentGradL2NormSquared ) ) << std::endl;
00590     }
00591     // currentGradInfNorm = 10e6;
00592     if( terminationCriterionFlag & ( GRADIENT_INF_NORM_ABSOLUTE | GRADIENT_INF_NORM_RELATIVE ) )
00593     {
00594         currentGradInfNorm = Linf( grad_array, pd.num_free_vertices() );  // get the Linf norm
00595         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- gradient Inf norm: "
00596                                          << " " << RPM( currentGradInfNorm ) << std::endl;
00597     }
00598 
00599     if( terminationCriterionFlag & VERTEX_MOVEMENT_RELATIVE )
00600     {
00601         maxSquaredInitialMovement = pd.get_max_vertex_movement_squared( initialVerticesMemento, err );MSQ_ERRRTN( err );
00602         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- max initial vertex movement: "
00603                                          << " " << RPM( maxSquaredInitialMovement ) << std::endl;
00604     }
00605 
00606     previousOFValue = currentOFValue;
00607     currentOFValue  = of_value;
00608     if( terminationCriterionFlag & OF_FLAGS )
00609     {
00610         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- OF Value: "
00611                                          << " " << RPM( of_value ) << " iterationCounter= " << iterationCounter
00612                                          << std::endl;
00613     }
00614     else if( grad_array )
00615     {
00616         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o OF Value: "
00617                                          << " " << RPM( of_value ) << " iterationCounter= "
00618                                          << iterationCounter
00619                                          //<< " terminationCriterionFlag= " <<
00620                                          // terminationCriterionFlag << " OF_FLAGS = " << OF_FLAGS
00621                                          << std::endl;
00622     }
00623 
00624     ++iterationCounter;
00625     if( timeStepFileType ) write_timestep( pd, grad_array, err );
00626 
00627     if( plotFile.is_open() )
00628         plotFile << iterationCounter << '\t' << mTimer.since_birth() << '\t' << of_value << '\t'
00629                  << std::sqrt( currentGradL2NormSquared ) << '\t' << currentGradInfNorm << '\t'
00630                  << ( maxSquaredMovement > 0.0 ? std::sqrt( maxSquaredMovement ) : 0.0 ) << '\t' << globalInvertedCount
00631                  << std::endl;
00632 }
00633 
00634 void TerminationCriterion::accumulate_outer( Mesh* mesh,
00635                                              MeshDomain* domain,
00636                                              OFEvaluator& of_eval,
00637                                              const Settings* settings,
00638                                              MsqError& err )
00639 {
00640     PatchData global_patch;
00641     if( settings ) global_patch.attach_settings( settings );
00642 
00643     // if we need to fill out the global patch data object.
00644     if( ( terminationCriterionFlag & ( GRAD_FLAGS | OF_FLAGS | VERTEX_MOVEMENT_RELATIVE ) ) || timeStepFileType )
00645     {
00646         global_patch.set_mesh( mesh );
00647         global_patch.set_domain( domain );
00648         global_patch.fill_global_patch( err );MSQ_ERRRTN( err );
00649     }
00650 
00651     accumulate_inner( global_patch, of_eval, err );MSQ_ERRRTN( err );
00652 }
00653 
00654 void TerminationCriterion::accumulate_patch( PatchData& pd, MsqError& err )
00655 {
00656     if( terminationCriterionFlag & MOVEMENT_FLAGS )
00657     {
00658         double patch_max_dist = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
00659         if( patch_max_dist > maxSquaredMovement ) maxSquaredMovement = patch_max_dist;
00660         pd.recreate_vertices_memento( previousVerticesMemento, err );MSQ_ERRRTN( err );
00661     }
00662 
00663     // if terminating on bounded vertex movement (a bounding box for the mesh)
00664     if( terminationCriterionFlag & BOUNDED_VERTEX_MOVEMENT )
00665     {
00666         const MsqVertex* vert = pd.get_vertex_array( err );
00667         int num_vert          = pd.num_free_vertices();
00668         int i                 = 0;
00669         // for each vertex
00670         for( i = 0; i < num_vert; ++i )
00671         {
00672             // if any of the coordinates are greater than eps
00673             if( ( vert[i][0] > boundedVertexMovementEps ) || ( vert[i][1] > boundedVertexMovementEps ) ||
00674                 ( vert[i][2] > boundedVertexMovementEps ) )
00675             {
00676                 ++vertexMovementExceedsBound;
00677             }
00678         }
00679     }
00680 
00681     if( ( terminationCriterionFlag | cullingMethodFlag ) & UNTANGLED_MESH )
00682     {
00683         size_t new_count = count_inverted( pd, err );
00684         // be careful here because size_t is unsigned
00685         globalInvertedCount += new_count;
00686         globalInvertedCount -= patchInvertedCount;
00687         patchInvertedCount = new_count;
00688         // if (innerOuterType==TYPE_OUTER)
00689         //  MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Patch Inverted: " << " " <<
00690         //  patchInvertedCount << " globalInvertedCount= " << globalInvertedCount << std::endl;
00691 
00692         MSQ_ERRRTN( err );
00693     }
00694 }
00695 
00696 /*!  This function evaluates the needed information and then evaluates
00697   the termination criteria.  If any of the selected criteria are satisfied,
00698   the function returns true.  Otherwise, the function returns false.
00699  */
00700 bool TerminationCriterion::terminate()
00701 {
00702     bool return_flag = false;
00703     // std::cout<<"\nInside terminate(pd,of,err):  flag = "<<terminationCriterionFlag << std::endl;
00704     int type = 0;
00705 
00706     // First check for an interrupt signal
00707     if( MsqInterrupt::interrupt() )
00708     {
00709         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- INTERRUPTED"
00710                                          << " " << std::endl;
00711         return true;
00712     }
00713 
00714     // if terminating on numbering of inner iterations
00715     if( NUMBER_OF_ITERATES & terminationCriterionFlag && iterationCounter >= iterationBound )
00716     {
00717         return_flag = true;
00718         type        = 1;
00719         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- Reached "
00720                                          << " " << iterationBound << " iterations." << std::endl;
00721     }
00722 
00723     if( CPU_TIME & terminationCriterionFlag && mTimer.since_birth() >= timeBound )
00724     {
00725         return_flag = true;
00726         type        = 2;
00727         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- Exceeded CPU time. "
00728                                          << " " << std::endl;
00729     }
00730 
00731     if( MOVEMENT_FLAGS & terminationCriterionFlag && maxSquaredMovement >= 0.0 )
00732     {
00733         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o Info -- Maximum vertex movement: "
00734                                          << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00735 
00736         if( VERTEX_MOVEMENT_ABSOLUTE & terminationCriterionFlag && maxSquaredMovement <= vertexMovementAbsoluteEps )
00737         {
00738             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_ABSOLUTE: "
00739                                              << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00740             return_flag = true;
00741             type        = 3;
00742         }
00743 
00744         if( VERTEX_MOVEMENT_RELATIVE & terminationCriterionFlag &&
00745             maxSquaredMovement <= vertexMovementRelativeEps * maxSquaredInitialMovement )
00746         {
00747             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_RELATIVE: "
00748                                              << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00749             return_flag = true;
00750             type        = 4;
00751         }
00752 
00753         if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & terminationCriterionFlag )
00754         {
00755             assert( vertexMovementAbsoluteAvgEdge > -1e-12 );  // make sure value actually got calculated
00756             if( maxSquaredMovement <= vertexMovementAbsoluteAvgEdge )
00757             {
00758                 MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- VERTEX_MOVEMENT_ABS_EDGE_LENGTH: "
00759                                                  << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
00760                 return_flag = true;
00761                 type        = 5;
00762             }
00763         }
00764     }
00765 
00766     if( GRADIENT_L2_NORM_ABSOLUTE & terminationCriterionFlag &&
00767         currentGradL2NormSquared <= gradL2NormAbsoluteEpsSquared )
00768     {
00769         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_L2_NORM_ABSOLUTE: "
00770                                          << " " << RPM( currentGradL2NormSquared ) << std::endl;
00771         return_flag = true;
00772         type        = 6;
00773     }
00774 
00775     if( GRADIENT_INF_NORM_ABSOLUTE & terminationCriterionFlag && currentGradInfNorm <= gradInfNormAbsoluteEps )
00776     {
00777         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_INF_NORM_ABSOLUTE: "
00778                                          << " " << RPM( currentGradInfNorm ) << std::endl;
00779         return_flag = true;
00780         type        = 7;
00781     }
00782 
00783     if( GRADIENT_L2_NORM_RELATIVE & terminationCriterionFlag &&
00784         currentGradL2NormSquared <= ( gradL2NormRelativeEpsSquared * initialGradL2NormSquared ) )
00785     {
00786         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_L2_NORM_RELATIVE: "
00787                                          << " " << RPM( currentGradL2NormSquared ) << std::endl;
00788         return_flag = true;
00789         type        = 8;
00790     }
00791 
00792     if( GRADIENT_INF_NORM_RELATIVE & terminationCriterionFlag &&
00793         currentGradInfNorm <= ( gradInfNormRelativeEps * initialGradInfNorm ) )
00794     {
00795         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- GRADIENT_INF_NORM_RELATIVE: "
00796                                          << " " << RPM( currentGradInfNorm ) << std::endl;
00797         return_flag = true;
00798         type        = 9;
00799     }
00800     // Quality Improvement and Successive Improvements are below.
00801     // The relative forms are only valid after the first iteration.
00802     if( ( QUALITY_IMPROVEMENT_ABSOLUTE & terminationCriterionFlag ) && currentOFValue <= qualityImprovementAbsoluteEps )
00803     {
00804         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- QUALITY_IMPROVEMENT_ABSOLUTE: "
00805                                          << " " << RPM( currentOFValue ) << std::endl;
00806         return_flag = true;
00807         type        = 10;
00808     }
00809 
00810     // only valid if an iteration has occurred, see above.
00811     if( iterationCounter > 0 )
00812     {
00813         if( SUCCESSIVE_IMPROVEMENTS_ABSOLUTE & terminationCriterionFlag &&
00814             ( previousOFValue - currentOFValue ) <= successiveImprovementsAbsoluteEps )
00815         {
00816             MSQ_DBGOUT_P0_ONLY( debugLevel )
00817                 << par_string() << "  o TermCrit -- SUCCESSIVE_IMPROVEMENTS_ABSOLUTE: previousOFValue= "
00818                 << " " << previousOFValue << " currentOFValue= " << RPM( currentOFValue )
00819                 << " successiveImprovementsAbsoluteEps= " << successiveImprovementsAbsoluteEps << std::endl;
00820             return_flag = true;
00821             type        = 11;
00822         }
00823         if( QUALITY_IMPROVEMENT_RELATIVE & terminationCriterionFlag &&
00824             ( currentOFValue - lowerOFBound ) <= qualityImprovementRelativeEps * ( initialOFValue - lowerOFBound ) )
00825         {
00826             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- QUALITY_IMPROVEMENT_RELATIVE: "
00827                                              << " " << std::endl;
00828             return_flag = true;
00829             type        = 12;
00830         }
00831         if( SUCCESSIVE_IMPROVEMENTS_RELATIVE & terminationCriterionFlag &&
00832             ( previousOFValue - currentOFValue ) <=
00833                 successiveImprovementsRelativeEps * ( initialOFValue - currentOFValue ) )
00834         {
00835             MSQ_DBGOUT_P0_ONLY( debugLevel )
00836                 << par_string() << "  o TermCrit -- SUCCESSIVE_IMPROVEMENTS_RELATIVE: previousOFValue= "
00837                 << " " << previousOFValue << " currentOFValue= " << RPM( currentOFValue )
00838                 << " successiveImprovementsRelativeEps= " << successiveImprovementsRelativeEps << std::endl;
00839             return_flag = true;
00840             type        = 13;
00841         }
00842     }
00843 
00844     if( BOUNDED_VERTEX_MOVEMENT & terminationCriterionFlag && vertexMovementExceedsBound )
00845     {
00846         return_flag = true;
00847         type        = 14;
00848         MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- "
00849                                          << " " << vertexMovementExceedsBound << " vertices out of bounds."
00850                                          << std::endl;
00851     }
00852 
00853     if( UNTANGLED_MESH & terminationCriterionFlag )
00854     {
00855         if( innerOuterType == TYPE_OUTER )
00856         {
00857             // MSQ_DBGOUT_P0_ONLY(debugLevel)
00858             MSQ_DBGOUT( debugLevel ) << par_string() << "  o Num Inverted: "
00859                                      << " " << globalInvertedCount << std::endl;
00860         }
00861         if( !globalInvertedCount )
00862         {
00863             MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << "  o TermCrit -- UNTANGLED_MESH: "
00864                                              << " " << std::endl;
00865             return_flag = true;
00866             type        = 15;
00867         }
00868     }
00869 
00870     // clear this value at the end of each iteration
00871     vertexMovementExceedsBound = 0;
00872     maxSquaredMovement         = -1.0;
00873 
00874     if( timeStepFileType == GNUPLOT && return_flag )
00875     {
00876         MsqError err;
00877         MeshWriter::write_gnuplot_overlay( iterationCounter, timeStepFileName.c_str(), err );
00878     }
00879 
00880     if( 0 && return_flag && MSQ_DBG( 2 ) )
00881         std::cout << "P[" << get_parallel_rank() << "] tmp TerminationCriterion::terminate: " << moniker
00882                   << " return_flag= " << return_flag << " type= " << type
00883                   << " terminationCriterionFlag= " << terminationCriterionFlag << " debugLevel= " << debugLevel
00884                   << std::endl;
00885 
00886     // if none of the criteria were satisfied
00887     return return_flag;
00888 }
00889 
00890 bool TerminationCriterion::criterion_is_set()
00891 {
00892     if( !terminationCriterionFlag )
00893         return false;
00894     else
00895         return true;
00896 }
00897 
00898 /*!This function checks the culling method criterion supplied to the object
00899   by the user.  If the user does not supply a culling method criterion,
00900   the default criterion is NONE, and in that case, no culling is performed.
00901   If the culling method criterion is satisfied, the interior vertices
00902   of the given patch are flagged as soft_fixed.  Otherwise, the soft_fixed
00903   flag is removed from each of the vertices in the patch (interior and
00904   boundary vertices).  Also, if the criterion was satisfied, then the
00905   function returns true.  Otherwise, the function returns false.
00906  */
00907 bool TerminationCriterion::cull_vertices( PatchData& pd, OFEvaluator& of_eval, MsqError& err )
00908 {
00909     // PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
00910 
00911     // cull_bool will be changed to true if the criterion is satisfied
00912     bool b, cull_bool = false;
00913     double prev_m, init_m;
00914     switch( cullingMethodFlag )
00915     {
00916             // if no culling is requested, always return false
00917         case NONE:
00918             return cull_bool;
00919             // if culling on quality improvement absolute
00920         case QUALITY_IMPROVEMENT_ABSOLUTE:
00921             // get objective function value
00922             b = of_eval.evaluate( pd, currentOFValue, err );
00923             if( MSQ_CHKERR( err ) ) return false;
00924             if( !b )
00925             {
00926                 MSQ_SETERR( err )( MsqError::INVALID_MESH );
00927                 return false;
00928             }
00929             // if the improvement was enough, cull
00930             if( currentOFValue <= cullingEps )
00931             {
00932                 cull_bool = true;
00933             }
00934             // PRINT_INFO("\ncurrentOFValue = %f, bool = %i\n",currentOFValue,cull_bool);
00935 
00936             break;
00937             // if culing on quality improvement relative
00938         case QUALITY_IMPROVEMENT_RELATIVE:
00939             // get objective function value
00940             b = of_eval.evaluate( pd, currentOFValue, err );
00941             if( MSQ_CHKERR( err ) ) return false;
00942             if( !b )
00943             {
00944                 MSQ_SETERR( err )( MsqError::INVALID_MESH );
00945                 return false;
00946             }
00947             // if the improvement was enough, cull
00948             if( ( currentOFValue - lowerOFBound ) <= ( cullingEps * ( initialOFValue - lowerOFBound ) ) )
00949             {
00950                 cull_bool = true;
00951             }
00952             break;
00953             // if culling on vertex movement absolute
00954         case VERTEX_MOVEMENT_ABSOLUTE:
00955         case VERTEX_MOVEMENT_ABS_EDGE_LENGTH:
00956             // if movement was enough, cull
00957             prev_m = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
00958             MSQ_ERRZERO( err );
00959             if( prev_m <= cullingEps * cullingEps )
00960             {
00961                 cull_bool = true;
00962             }
00963 
00964             break;
00965             // if culling on vertex movement relative
00966         case VERTEX_MOVEMENT_RELATIVE:
00967             // if movement was small enough, cull
00968             prev_m = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
00969             MSQ_ERRZERO( err );
00970             init_m = pd.get_max_vertex_movement_squared( initialVerticesMemento, err );
00971             MSQ_ERRZERO( err );
00972             if( prev_m <= ( cullingEps * cullingEps * init_m ) )
00973             {
00974                 cull_bool = true;
00975             }
00976             break;
00977         case UNTANGLED_MESH:
00978             if( !patchInvertedCount ) cull_bool = true;
00979             break;
00980         default:
00981             MSQ_SETERR( err )
00982             ( "Requested culling method not yet implemented.", MsqError::NOT_IMPLEMENTED );
00983             return false;
00984     };
00985     // Now actually have patch data cull vertices
00986     if( cull_bool )
00987     {
00988         pd.set_free_vertices_soft_fixed( err );
00989         MSQ_ERRZERO( err );
00990     }
00991     else
00992     {
00993         pd.set_all_vertices_soft_free( err );
00994         MSQ_ERRZERO( err );
00995     }
00996     return cull_bool;
00997 }
00998 
00999 /*!This function is activated when cullingGlobalPatch is true.  It supplies
01000   cull_vertices with a single vertex-based patch at a time.  If the patch
01001   satisfies the culling criterion, it's free vertices are then soft-fixed.
01002  */
01003 bool TerminationCriterion::cull_vertices_global( PatchData& global_patch,
01004                                                  Mesh* mesh,
01005                                                  MeshDomain* domain,
01006                                                  const Settings* settings,
01007                                                  OFEvaluator& of_eval,
01008                                                  MsqError& err )
01009 {
01010     if( !cullingGlobalPatch ) return false;
01011 
01012     // PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
01013 
01014     // cull_bool will be changed to true if the criterion is satisfied
01015     bool cull_bool = false;
01016 
01017     std::vector< Mesh::VertexHandle > mesh_vertices;
01018     // std::vector<Mesh::VertexHandle> patch_vertices;
01019     // std::vector<Mesh::ElementHandle> patch_elements;
01020     // std::vector<Mesh::VertexHandle> fixed_vertices;
01021     // std::vector<Mesh::VertexHandle> free_vertices;
01022 
01023     // FIXME, verify global_patch is a global patch... how, is this right?
01024     mesh->get_all_vertices( mesh_vertices, err );
01025     size_t mesh_num_nodes         = mesh_vertices.size();
01026     size_t global_patch_num_nodes = global_patch.num_nodes();
01027     if( 0 )
01028         std::cout << "tmp srk mesh_num_nodes= " << mesh_num_nodes
01029                   << " global_patch_num_nodes= " << global_patch_num_nodes << std::endl;
01030     if( mesh_num_nodes != global_patch_num_nodes )
01031     {
01032         std::cout << "tmp srk cull_vertices_global found non global patch" << std::endl;
01033         exit( 123 );
01034         return false;
01035     }
01036     PatchData patch;
01037     patch.set_mesh( (Mesh*)mesh );
01038     patch.set_domain( domain );
01039     patch.attach_settings( settings );
01040 
01041     // const MsqVertex* global_patch_vertex_array = global_patch.get_vertex_array( err );
01042     Mesh::VertexHandle* global_patch_vertex_handles = global_patch.get_vertex_handles_array();
01043 
01044     int num_culled = 0;
01045     for( unsigned iv = 0; iv < global_patch_num_nodes; iv++ )
01046     {
01047         // form a patch for this vertex; if it is culled, set it to be soft fixed
01048         Mesh::VertexHandle vert = global_patch_vertex_handles[iv];
01049         std::vector< Mesh::ElementHandle > elements;
01050         std::vector< size_t > offsets;
01051         mesh->vertices_get_attached_elements( &vert, 1, elements, offsets, err );
01052 
01053         std::set< Mesh::VertexHandle > patch_free_vertices_set;
01054 
01055         for( unsigned ie = 0; ie < elements.size(); ie++ )
01056         {
01057             std::vector< Mesh::VertexHandle > vert_handles;
01058             std::vector< size_t > v_offsets;
01059             mesh->elements_get_attached_vertices( &elements[ie], 1, vert_handles, v_offsets, err );
01060             for( unsigned jv = 0; jv < vert_handles.size(); jv++ )
01061             {
01062                 unsigned char bt;
01063                 mesh->vertex_get_byte( vert_handles[jv], &bt, err );
01064                 MsqVertex v;
01065                 v.set_flags( bt );
01066                 if( v.is_free_vertex() ) patch_free_vertices_set.insert( vert_handles[jv] );
01067             }
01068         }
01069 
01070         std::vector< Mesh::VertexHandle > patch_free_vertices_vector( patch_free_vertices_set.begin(),
01071                                                                       patch_free_vertices_set.end() );
01072         // std::vector<unsigned char> byte_vector(patch_vertices_vector.size());
01073         // mesh->vertices_get_byte(&vert_handles[0], &byte_vector[0], vert_handles.size(), err);
01074 
01075         patch.set_mesh_entities( elements, patch_free_vertices_vector, err );
01076         if( cull_vertices( patch, of_eval, err ) )
01077         {
01078             // std::cout << "tmp srk cull_vertices_global found culled patch" << std::endl;
01079             Mesh::VertexHandle* patch_vertex_handles = patch.get_vertex_handles_array();
01080             const MsqVertex* patch_vertex_array      = patch.get_vertex_array( err );
01081             for( unsigned jv = 0; jv < patch.num_nodes(); jv++ )
01082             {
01083                 if( patch_vertex_handles[jv] == global_patch_vertex_handles[iv] )
01084                 {
01085                     if( patch_vertex_array[jv].is_flag_set( MsqVertex::MSQ_CULLED ) )
01086                     {
01087                         global_patch.set_vertex_culled( iv );
01088                         ++num_culled;
01089                         cull_bool = true;
01090                         // std::cout << "tmp srk cull_vertices_global found culled vertex" <<
01091                         // std::endl;
01092                     }
01093                 }
01094             }
01095         }
01096     }
01097     if( 0 )
01098         std::cout << "tmp srk cull_vertices_global found " << num_culled << " culled vertices out of "
01099                   << global_patch_num_nodes << std::endl;
01100 
01101     return cull_bool;
01102 }
01103 
01104 size_t TerminationCriterion::count_inverted( PatchData& pd, MsqError& err )
01105 {
01106     size_t num_elem = pd.num_elements();
01107     size_t count    = 0;
01108     int inverted, samples;
01109     for( size_t i = 0; i < num_elem; i++ )
01110     {
01111         pd.element_by_index( i ).check_element_orientation( pd, inverted, samples, err );
01112         if( inverted ) ++count;
01113     }
01114     return count;
01115 }
01116 
01117 /*!
01118   Currently this only deletes the memento of the vertex positions and the
01119   mGrad vector if neccessary.
01120   When culling, we remove the soft fixed flags from all of the vertices.
01121  */
01122 void TerminationCriterion::cleanup( Mesh*, MeshDomain*, MsqError& )
01123 {
01124     delete previousVerticesMemento;
01125     delete initialVerticesMemento;
01126     previousVerticesMemento = 0;
01127     initialVerticesMemento  = 0;
01128 }
01129 
01130 void TerminationCriterion::write_timestep( PatchData& pd, const Vector3D* gradient, MsqError& err )
01131 {
01132     std::ostringstream str;
01133     if( timeStepFileType == VTK )
01134     {
01135         str << timeStepFileName << '_' << iterationCounter << ".vtk";
01136         MeshWriter::write_vtk( pd, str.str().c_str(), err, gradient );
01137     }
01138     else if( timeStepFileType == GNUPLOT )
01139     {
01140         str << timeStepFileName << '.' << iterationCounter;
01141         MeshWriter::write_gnuplot( pd.get_mesh(), str.str().c_str(), err );
01142     }
01143 }
01144 
01145 void TerminationCriterion::write_iterations( const char* filename, MsqError& err )
01146 {
01147     if( filename )
01148     {
01149         plotFile.open( filename, std::ios::trunc );
01150         if( !plotFile ) MSQ_SETERR( err )
01151         ( MsqError::FILE_ACCESS, "Failed to open plot data file: '%s'", filename );
01152     }
01153     else
01154     {
01155         plotFile.close();
01156     }
01157 }
01158 
01159 void TerminationCriterion::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings*, MsqError& err )
01160 {
01161     if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & ( terminationCriterionFlag | cullingMethodFlag ) )
01162     {
01163         Mesh* mesh = mesh_and_domain->get_mesh();
01164         MeshUtil tool( mesh );
01165         SimpleStats stats;
01166         tool.edge_length_distribution( stats, err );MSQ_ERRRTN( err );
01167         double limit = vertexMovementAvgBeta * ( stats.average() - stats.standard_deviation() );
01168         // we actually calculate the square of the length
01169         vertexMovementAbsoluteAvgEdge = limit * limit;
01170         if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & cullingMethodFlag ) cullingEps = limit;
01171     }
01172 }
01173 
01174 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines