Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : :
26 : : ***************************************************************** */
27 : : // -*- Mode : c++; tab-width: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2
28 : : // -*-
29 : : //
30 : : // DESCRIPTION:
31 : : // ============
32 : : /*! \file TerminationCriterion.cpp
33 : :
34 : : \brief Member functions of the MBMesquite::TerminationCriterion class
35 : :
36 : : \author Michael Brewer
37 : : \author Thomas Leurent
38 : : \date Feb. 14, 2003
39 : : */
40 : :
41 : : #include "TerminationCriterion.hpp"
42 : : #include "MsqVertex.hpp"
43 : : #include "MsqInterrupt.hpp"
44 : : #include "OFEvaluator.hpp"
45 : : #include "MsqError.hpp"
46 : : #include "MsqDebug.hpp"
47 : : #include "PatchData.hpp"
48 : : #include "MeshWriter.hpp"
49 : : #include "MeshUtil.hpp"
50 : : #include "SimpleStats.hpp"
51 : :
52 : : #include <sstream>
53 : : #include <set>
54 : :
55 : : namespace MBMesquite
56 : : {
57 : :
58 : : extern int get_parallel_rank();
59 : : extern int get_parallel_size();
60 : : extern double reduce_parallel_max( double value );
61 : :
62 : : #define MSQ_DBGOUT_P0_ONLY( flag ) \
63 : : if( !get_parallel_rank() ) MSQ_DBGOUT( flag )
64 : :
65 : : #define RPM( val ) val
66 : :
67 : : // this causes race conditions - don't use it
68 : : #define RPM1( val ) reduce_parallel_max( val )
69 : :
70 : : /*! \enum TCType defines the termination criterion */
71 : : enum TCType
72 : : {
73 : : NONE = 0,
74 : : //! checks the gradient \f$\nabla f \f$ of objective function
75 : : //! \f$f : I\!\!R^{3N} \rightarrow I\!\!R \f$ against a double \f$d\f$
76 : : //! and stops when \f$\sqrt{\sum_{i=1}^{3N}\nabla f_i^2}<d\f$
77 : : GRADIENT_L2_NORM_ABSOLUTE = 1 << 0,
78 : : //! checks the gradient \f$\nabla f \f$ of objective function
79 : : //! \f$f : I\!\!R^{3N} \rightarrow I\!\!R \f$ against a double \f$d\f$
80 : : //! and stops when \f$ \max_{i=1}^{3N} \nabla f_i < d \f$
81 : : GRADIENT_INF_NORM_ABSOLUTE = 1 << 1,
82 : : //! terminates on the j_th iteration when
83 : : //! \f$\sqrt{\sum_{i=1}^{3N}\nabla f_{i,j}^2}<d\sqrt{\sum_{i=1}^{3N}\nabla f_{i,0}^2}\f$
84 : : //! That is, terminates when the norm of the gradient is small
85 : : //! than some scaling factor times the norm of the original gradient.
86 : : GRADIENT_L2_NORM_RELATIVE = 1 << 2,
87 : : //! terminates on the j_th iteration when
88 : : //! \f$\max_{i=1 \cdots 3N}\nabla f_{i,j}<d \max_{i=1 \cdots 3N}\nabla f_{i,0}\f$
89 : : //! That is, terminates when the norm of the gradient is small
90 : : //! than some scaling factor times the norm of the original gradient.
91 : : //! (Using the infinity norm.)
92 : : GRADIENT_INF_NORM_RELATIVE = 1 << 3,
93 : : //! Not yet implemented.
94 : : KKT = 1 << 4,
95 : : //! Terminates when the objective function value is smaller than
96 : : //! the given scalar value.
97 : : QUALITY_IMPROVEMENT_ABSOLUTE = 1 << 5,
98 : : //! Terminates when the objective function value is smaller than
99 : : //! the given scalar value times the original objective function
100 : : //! value.
101 : : QUALITY_IMPROVEMENT_RELATIVE = 1 << 6,
102 : : //! Terminates when the number of iterations exceeds a given integer.
103 : : NUMBER_OF_ITERATES = 1 << 7,
104 : : //! Terminates when the algorithm exceeds an allotted time limit
105 : : //! (given in seconds).
106 : : CPU_TIME = 1 << 8,
107 : : //! Terminates when a the maximum distance moved by any vertex
108 : : //! during the previous iteration is below the given value.
109 : : VERTEX_MOVEMENT_ABSOLUTE = 1 << 9,
110 : : //! Terminates when a the maximum distance moved by any vertex
111 : : //! during the previous iteration is below a value calculated
112 : : //! from the edge lengths of the initial mesh.
113 : : VERTEX_MOVEMENT_ABS_EDGE_LENGTH = 1 << 10,
114 : : //! Terminates when a the maximum distance moved by any vertex
115 : : //! during the previous iteration is below the given value
116 : : //! times the maximum distance moved by any vertex over the
117 : : //! entire course of the optimization.
118 : : VERTEX_MOVEMENT_RELATIVE = 1 << 11,
119 : : //! Terminates when the decrease in the objective function value since
120 : : //! the previous iteration is below the given value.
121 : : SUCCESSIVE_IMPROVEMENTS_ABSOLUTE = 1 << 12,
122 : : //! Terminates when the decrease in the objective function value since
123 : : //! the previous iteration is below the given value times the
124 : : //! decrease in the objective function value since the beginning
125 : : //! of this optimization process.
126 : : SUCCESSIVE_IMPROVEMENTS_RELATIVE = 1 << 13,
127 : : //! Terminates when any vertex leaves the bounding box, defined
128 : : //! by the given value, d. That is, when the absolute value of
129 : : //! a single coordinate of vertex's position exceeds d.
130 : : BOUNDED_VERTEX_MOVEMENT = 1 << 14,
131 : : //! Terminate when no elements are inverted
132 : : UNTANGLED_MESH = 1 << 15
133 : : };
134 : :
135 : : const unsigned long GRAD_FLAGS =
136 : : GRADIENT_L2_NORM_ABSOLUTE | GRADIENT_INF_NORM_ABSOLUTE | GRADIENT_L2_NORM_RELATIVE | GRADIENT_INF_NORM_RELATIVE;
137 : : const unsigned long OF_FLAGS = QUALITY_IMPROVEMENT_ABSOLUTE | QUALITY_IMPROVEMENT_RELATIVE |
138 : : SUCCESSIVE_IMPROVEMENTS_ABSOLUTE | SUCCESSIVE_IMPROVEMENTS_RELATIVE;
139 : :
140 : : const unsigned long MOVEMENT_FLAGS =
141 : : VERTEX_MOVEMENT_ABSOLUTE | VERTEX_MOVEMENT_ABS_EDGE_LENGTH | VERTEX_MOVEMENT_RELATIVE;
142 : :
143 : : /*!Constructor initializes all of the data members which are not
144 : : necessarily automatically initialized in their constructors.*/
145 : 569 : TerminationCriterion::TerminationCriterion( std::string name, InnerOuterType pinnerOuterType )
146 : : : mGrad( 8 ), initialVerticesMemento( 0 ), previousVerticesMemento( 0 ), debugLevel( 2 ),
147 [ + - ][ + - ]: 569 : timeStepFileType( NOTYPE ), moniker( name ), innerOuterType( pinnerOuterType )
[ + - ][ + - ]
[ + - ]
148 : : {
149 : 569 : terminationCriterionFlag = NONE;
150 : 569 : cullingMethodFlag = NONE;
151 : 569 : cullingEps = 0.0;
152 : 569 : cullingGlobalPatch = false;
153 : 569 : initialOFValue = 0.0;
154 : 569 : previousOFValue = 0.0;
155 : 569 : currentOFValue = 0.0;
156 : 569 : lowerOFBound = 0.0;
157 : 569 : initialGradL2NormSquared = 0.0;
158 : 569 : initialGradInfNorm = 0.0;
159 : : // initial size of the gradient array
160 : 569 : gradL2NormAbsoluteEpsSquared = 0.0;
161 : 569 : gradL2NormRelativeEpsSquared = 0.0;
162 : 569 : gradInfNormAbsoluteEps = 0.0;
163 : 569 : gradInfNormRelativeEps = 0.0;
164 : 569 : qualityImprovementAbsoluteEps = 0.0;
165 : 569 : qualityImprovementRelativeEps = 0.0;
166 : 569 : iterationBound = 0;
167 : 569 : iterationCounter = 0;
168 : 569 : timeBound = 0.0;
169 : 569 : vertexMovementAbsoluteEps = 0.0;
170 : 569 : vertexMovementRelativeEps = 0.0;
171 : 569 : vertexMovementAvgBeta = -1;
172 : 569 : vertexMovementAbsoluteAvgEdge = -1;
173 : 569 : successiveImprovementsAbsoluteEps = 0.0;
174 : 569 : successiveImprovementsRelativeEps = 0.0;
175 : 569 : boundedVertexMovementEps = 0.0;
176 : 569 : globalInvertedCount = 0;
177 : 569 : patchInvertedCount = 0;
178 : 569 : }
179 : :
180 : 0 : std::string TerminationCriterion::par_string()
181 : : {
182 [ # # ]: 0 : if( get_parallel_size() )
183 : : {
184 [ # # ]: 0 : std::ostringstream str;
185 [ # # ][ # # ]: 0 : str << "P[" << get_parallel_rank() << "] " + moniker + " ";
[ # # ][ # # ]
[ # # ][ # # ]
186 [ # # ]: 0 : return str.str();
187 : : }
188 : 0 : return moniker + " ";
189 : : }
190 : :
191 : 59 : void TerminationCriterion::add_absolute_gradient_L2_norm( double eps )
192 : : {
193 : 59 : terminationCriterionFlag |= GRADIENT_L2_NORM_ABSOLUTE;
194 : 59 : gradL2NormAbsoluteEpsSquared = eps * eps;
195 : 59 : }
196 : :
197 : 0 : void TerminationCriterion::add_absolute_gradient_inf_norm( double eps )
198 : : {
199 : 0 : terminationCriterionFlag |= GRADIENT_INF_NORM_ABSOLUTE;
200 : 0 : gradInfNormAbsoluteEps = eps;
201 : 0 : }
202 : :
203 : 0 : void TerminationCriterion::add_relative_gradient_L2_norm( double eps )
204 : : {
205 : 0 : terminationCriterionFlag |= GRADIENT_L2_NORM_RELATIVE;
206 : 0 : gradL2NormRelativeEpsSquared = eps * eps;
207 : 0 : }
208 : :
209 : 0 : void TerminationCriterion::add_relative_gradient_inf_norm( double eps )
210 : : {
211 : 0 : terminationCriterionFlag |= GRADIENT_INF_NORM_RELATIVE;
212 : 0 : gradInfNormRelativeEps = eps;
213 : 0 : }
214 : :
215 : 3 : void TerminationCriterion::add_absolute_quality_improvement( double eps )
216 : : {
217 : 3 : terminationCriterionFlag |= QUALITY_IMPROVEMENT_ABSOLUTE;
218 : 3 : qualityImprovementAbsoluteEps = eps;
219 : 3 : }
220 : :
221 : 2 : void TerminationCriterion::add_relative_quality_improvement( double eps )
222 : : {
223 : 2 : terminationCriterionFlag |= QUALITY_IMPROVEMENT_RELATIVE;
224 : 2 : qualityImprovementRelativeEps = eps;
225 : 2 : }
226 : :
227 : 79 : void TerminationCriterion::add_absolute_vertex_movement( double eps )
228 : : {
229 : 79 : terminationCriterionFlag |= VERTEX_MOVEMENT_ABSOLUTE;
230 : : // we actually compare squared movement to squared epsilon
231 : 79 : vertexMovementAbsoluteEps = ( eps * eps );
232 : 79 : }
233 : :
234 : 3 : void TerminationCriterion::add_absolute_vertex_movement_edge_length( double beta )
235 : : {
236 : 3 : terminationCriterionFlag |= VERTEX_MOVEMENT_ABS_EDGE_LENGTH;
237 : 3 : vertexMovementAvgBeta = beta;
238 : 3 : }
239 : :
240 : 0 : void TerminationCriterion::add_relative_vertex_movement( double eps )
241 : : {
242 : 0 : terminationCriterionFlag |= VERTEX_MOVEMENT_RELATIVE;
243 : : // we actually compare squared movement to squared epsilon
244 : 0 : vertexMovementRelativeEps = ( eps * eps );
245 : 0 : }
246 : :
247 : 2 : void TerminationCriterion::add_absolute_successive_improvement( double eps )
248 : : {
249 : 2 : terminationCriterionFlag |= SUCCESSIVE_IMPROVEMENTS_ABSOLUTE;
250 : 2 : successiveImprovementsAbsoluteEps = eps;
251 : 2 : }
252 : :
253 : 2 : void TerminationCriterion::add_relative_successive_improvement( double eps )
254 : : {
255 : 2 : terminationCriterionFlag |= SUCCESSIVE_IMPROVEMENTS_RELATIVE;
256 : 2 : successiveImprovementsRelativeEps = eps;
257 : 2 : }
258 : :
259 : 4 : void TerminationCriterion::add_cpu_time( double seconds )
260 : : {
261 : 4 : terminationCriterionFlag |= CPU_TIME;
262 : 4 : timeBound = seconds;
263 : 4 : }
264 : :
265 : 347 : void TerminationCriterion::add_iteration_limit( unsigned int max_iterations )
266 : : {
267 : 347 : terminationCriterionFlag |= NUMBER_OF_ITERATES;
268 : 347 : iterationBound = max_iterations;
269 : 347 : }
270 : :
271 : 0 : void TerminationCriterion::add_bounded_vertex_movement( double eps )
272 : : {
273 : 0 : terminationCriterionFlag |= BOUNDED_VERTEX_MOVEMENT;
274 : 0 : boundedVertexMovementEps = eps;
275 : 0 : }
276 : :
277 : 20 : void TerminationCriterion::add_untangled_mesh()
278 : : {
279 : 20 : terminationCriterionFlag |= UNTANGLED_MESH;
280 : 20 : }
281 : :
282 : 0 : void TerminationCriterion::remove_all_criteria()
283 : : {
284 : 0 : terminationCriterionFlag = 0;
285 : 0 : }
286 : :
287 : 0 : void TerminationCriterion::cull_on_absolute_quality_improvement( double limit )
288 : : {
289 : 0 : cullingMethodFlag = QUALITY_IMPROVEMENT_ABSOLUTE;
290 : 0 : cullingEps = limit;
291 : 0 : }
292 : 0 : void TerminationCriterion::cull_on_relative_quality_improvement( double limit )
293 : : {
294 : 0 : cullingMethodFlag = QUALITY_IMPROVEMENT_RELATIVE;
295 : 0 : cullingEps = limit;
296 : 0 : }
297 : 20 : void TerminationCriterion::cull_on_absolute_vertex_movement( double limit )
298 : : {
299 : 20 : cullingMethodFlag = VERTEX_MOVEMENT_ABSOLUTE;
300 : 20 : cullingEps = limit;
301 : 20 : }
302 : 0 : void TerminationCriterion::cull_on_relative_vertex_movement( double limit )
303 : : {
304 : 0 : cullingMethodFlag = VERTEX_MOVEMENT_RELATIVE;
305 : 0 : cullingEps = limit;
306 : 0 : }
307 : 7 : void TerminationCriterion::cull_on_absolute_vertex_movement_edge_length( double limit )
308 : : {
309 : 7 : cullingMethodFlag = VERTEX_MOVEMENT_ABS_EDGE_LENGTH;
310 : 7 : vertexMovementAvgBeta = limit;
311 : 7 : }
312 : 0 : void TerminationCriterion::cull_on_absolute_successive_improvement( double limit )
313 : : {
314 : 0 : cullingMethodFlag = SUCCESSIVE_IMPROVEMENTS_ABSOLUTE;
315 : 0 : cullingEps = limit;
316 : 0 : }
317 : 0 : void TerminationCriterion::cull_on_relative_successive_improvement( double limit )
318 : : {
319 : 0 : cullingMethodFlag = SUCCESSIVE_IMPROVEMENTS_RELATIVE;
320 : 0 : cullingEps = limit;
321 : 0 : }
322 : 0 : void TerminationCriterion::cull_untangled_mesh()
323 : : {
324 : 0 : cullingMethodFlag = UNTANGLED_MESH;
325 : 0 : }
326 : :
327 : 0 : void TerminationCriterion::remove_culling()
328 : : {
329 : 0 : cullingMethodFlag = NONE;
330 : 0 : }
331 : :
332 : 0 : void TerminationCriterion::cull_for_global_patch( bool val )
333 : : {
334 : 0 : cullingGlobalPatch = val;
335 : 0 : }
336 : :
337 : : /*!This version of reset is called using a MeshSet, which implies
338 : : it is only called when this criterion is used as the 'outer' termination
339 : : criterion.
340 : : */
341 : 128 : void TerminationCriterion::reset_outer( Mesh* mesh, MeshDomain* domain, OFEvaluator& obj_eval, const Settings* settings,
342 : : MsqError& err )
343 : : {
344 : 128 : const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
345 [ + - ]: 128 : PatchData global_patch;
346 [ + - ][ + - ]: 128 : if( settings ) global_patch.attach_settings( settings );
347 : :
348 : : // if we need to fill out the global patch data object.
349 [ + + ][ - + ]: 128 : if( ( totalFlag & ( GRAD_FLAGS | OF_FLAGS | VERTEX_MOVEMENT_RELATIVE | UNTANGLED_MESH ) ) || timeStepFileType )
350 : : {
351 [ + - ]: 21 : global_patch.set_mesh( mesh );
352 [ + - ]: 21 : global_patch.set_domain( domain );
353 [ + - ][ + - ]: 21 : global_patch.fill_global_patch( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
354 : : }
355 : :
356 : : // now call the other reset
357 [ + - ][ + - ]: 128 : reset_inner( global_patch, obj_eval, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
[ + - ]
358 : : }
359 : :
360 : : /*!Reset function using using a PatchData object. This function is
361 : : called for the inner-stopping criterion directly from the
362 : : loop over mesh function in VertexMover. For outer criterion,
363 : : it is called from the reset function which takes a MeshSet object.
364 : : This function prepares the object to be used by setting the initial
365 : : values of some of the data members. As examples, if needed, it resets
366 : : the cpu timer to zero, the iteration counter to zero, and the
367 : : initial and previous objective function values to the current
368 : : objective function value for this patch.
369 : : The return value for this function is similar to that of terminate().
370 : : The function returns false if the checked criteria have not been
371 : : satisfied, and true if they have been. reset() only checks the
372 : : GRADIENT_INF_NORM_ABSOLUTE, GRADIENT_L2_NORM_ABSOLUTE, and the
373 : : QUALITY_IMPROVEMENT_ABSOLUTE criteria. Checking these criteria
374 : : allows the QualityImprover to skip the entire optimization if
375 : : the initial mesh satisfies the appropriate conditions.
376 : : */
377 : 830076 : void TerminationCriterion::reset_inner( PatchData& pd, OFEvaluator& obj_eval, MsqError& err )
378 : : {
379 : 830076 : const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
380 : :
381 : : // clear flag for BOUNDED_VERTEX_MOVEMENT
382 : 830076 : vertexMovementExceedsBound = 0;
383 : :
384 : : // Use -1 to denote that this isn't initialized yet.
385 : : // As all valid values must be >= 0.0, a negative
386 : : // value indicates that it is uninitialized and is
387 : : // always less than any valid value.
388 : 830076 : maxSquaredMovement = -1;
389 : :
390 : : // Clear the iteration count.
391 : 830076 : iterationCounter = 0;
392 : :
393 : : // reset the inner timer if needed
394 [ + + ]: 830076 : if( totalFlag & CPU_TIME ) { mTimer.reset(); }
395 : :
396 : : // GRADIENT
397 : 830076 : currentGradInfNorm = initialGradInfNorm = 0.0;
398 : 830076 : currentGradL2NormSquared = initialGradL2NormSquared = 0.0;
399 [ + + ]: 830076 : if( totalFlag & GRAD_FLAGS )
400 : : {
401 [ - + ]: 3 : if( !obj_eval.have_objective_function() )
402 : : {
403 : : MSQ_SETERR( err )
404 : : ( "Error termination criteria set which uses objective "
405 : : "functions, but no objective function is available.",
406 [ # # ]: 0 : MsqError::INVALID_STATE );
407 : 0 : return;
408 : : }
409 : 3 : int num_vertices = pd.num_free_vertices();
410 : 3 : mGrad.resize( num_vertices );
411 : :
412 : : // get gradient and make sure it is valid
413 [ - + ][ # # ]: 3 : bool b = obj_eval.evaluate( pd, currentOFValue, mGrad, err );MSQ_ERRRTN( err );
[ - + ]
414 [ - + ]: 3 : if( !b )
415 : : {
416 : : MSQ_SETERR( err )
417 [ # # ]: 0 : ( "Initial patch is invalid for gradient computation.", MsqError::INVALID_STATE );
418 : 0 : return;
419 : : }
420 : :
421 : : // get the gradient norms
422 [ - + ]: 3 : if( totalFlag & ( GRADIENT_INF_NORM_ABSOLUTE | GRADIENT_INF_NORM_RELATIVE ) )
423 : : {
424 : 0 : currentGradInfNorm = initialGradInfNorm = Linf( mGrad );
425 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o Initial gradient Inf norm: "
426 : : << " " << RPM( initialGradInfNorm ) << std::endl;
427 : : }
428 : :
429 [ + - ]: 3 : if( totalFlag & ( GRADIENT_L2_NORM_ABSOLUTE | GRADIENT_L2_NORM_RELATIVE ) )
430 : : {
431 : 3 : currentGradL2NormSquared = initialGradL2NormSquared = length_squared( mGrad );
432 : 3 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o Initial gradient L2 norm: "
433 : : << " " << RPM( std::sqrt( initialGradL2NormSquared ) ) << std::endl;
434 : : }
435 : :
436 : : // the OFvalue comes for free, so save it
437 : 3 : previousOFValue = currentOFValue;
438 : 3 : initialOFValue = currentOFValue;
439 : : }
440 : : // find the initial objective function value if needed and not already
441 : : // computed. If we needed the gradient, we have the OF value for free.
442 : : // Also, if possible, get initial OF value if writing plot file. Solvers
443 : : // often supply the OF value for subsequent iterations so by calculating
444 : : // the initial value we can generate OF value plots.
445 [ + + - + ]: 1660141 : else if( ( totalFlag & OF_FLAGS ) ||
[ + + ]
446 [ # # ][ # # ]: 830068 : ( plotFile.is_open() && pd.num_free_vertices() && obj_eval.have_objective_function() ) )
447 : : {
448 : : // ensure the obj_ptr is not null
449 [ - + ]: 5 : if( !obj_eval.have_objective_function() )
450 : : {
451 : : MSQ_SETERR( err )
452 : : ( "Error termination criteria set which uses objective "
453 : : "functions, but no objective function is available.",
454 [ # # ]: 0 : MsqError::INVALID_STATE );
455 : 0 : return;
456 : : }
457 : :
458 [ - + ][ # # ]: 5 : bool b = obj_eval.evaluate( pd, currentOFValue, err );MSQ_ERRRTN( err );
[ - + ]
459 [ - + ]: 5 : if( !b )
460 : : {
461 : : MSQ_SETERR( err )
462 [ # # ]: 0 : ( "Initial patch is invalid for evaluation.", MsqError::INVALID_STATE );
463 : 0 : return;
464 : : }
465 : : // std::cout<<"\nReseting initial of value = "<<initialOFValue;
466 : 5 : previousOFValue = currentOFValue;
467 : 5 : initialOFValue = currentOFValue;
468 : : }
469 : :
470 [ + + ]: 830076 : if( totalFlag & ( GRAD_FLAGS | OF_FLAGS ) )
471 : 8 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o Initial OF value: "
472 : : << " " << RPM( initialOFValue ) << std::endl;
473 : :
474 : : // Store current vertex locations now, because we'll
475 : : // need them later to compare the current movement with.
476 [ - + ]: 830076 : if( totalFlag & VERTEX_MOVEMENT_RELATIVE )
477 : : {
478 [ # # ]: 0 : if( initialVerticesMemento ) { pd.recreate_vertices_memento( initialVerticesMemento, err ); }
479 : : else
480 : : {
481 : 0 : initialVerticesMemento = pd.create_vertices_memento( err );
482 : : }
483 [ # # ][ # # ]: 0 : MSQ_ERRRTN( err );
[ # # ]
484 : 0 : maxSquaredInitialMovement = DBL_MAX;
485 : : }
486 : : else
487 : : {
488 : 830076 : maxSquaredInitialMovement = 0;
489 : : }
490 : :
491 [ + + ]: 830076 : if( terminationCriterionFlag & UNTANGLED_MESH )
492 : : {
493 : 20 : globalInvertedCount = count_inverted( pd, err );
494 : : // if (innerOuterType==TYPE_OUTER) MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << " o
495 : : // Num Inverted: " << " " << globalInvertedCount << std::endl;
496 [ - + ][ # # ]: 20 : patchInvertedCount = 0;MSQ_ERRRTN( err );
[ - + ]
497 : : }
498 : :
499 [ + + ]: 830076 : if( timeStepFileType )
500 : : {
501 : : // If didn't already calculate gradient abive, calculate it now.
502 [ + - ]: 1 : if( !( totalFlag & GRAD_FLAGS ) )
503 : : {
504 : 1 : mGrad.resize( pd.num_free_vertices() );
505 : 1 : obj_eval.evaluate( pd, currentOFValue, mGrad, err );
506 : 1 : err.clear();
507 : : }
508 [ - + ]: 1 : write_timestep( pd, mGrad.empty() ? 0 : arrptr( mGrad ), err );
509 : : }
510 : :
511 [ - + ]: 830076 : if( plotFile.is_open() )
512 : : {
513 : : // two newlines so GNU plot knows that we are starting a new data set
514 : 0 : plotFile << std::endl << std::endl;
515 : : // write column headings as comment in data file
516 : 0 : plotFile << "#Iter\tCPU\tObjFunc\tGradL2\tGradInf\tMovement\tInverted" << std::endl;
517 : : // write initial values
518 : 0 : plotFile << 0 << '\t' << mTimer.since_birth() << '\t' << initialOFValue << '\t'
519 : 0 : << std::sqrt( currentGradL2NormSquared ) << '\t' << currentGradInfNorm << '\t' << 0.0 << '\t'
520 : 830076 : << globalInvertedCount << std::endl;
521 : : }
522 : : }
523 : :
524 : 1659896 : void TerminationCriterion::reset_patch( PatchData& pd, MsqError& err )
525 : : {
526 : 1659896 : const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
527 [ + + ]: 1659896 : if( totalFlag & MOVEMENT_FLAGS )
528 : : {
529 [ + + ]: 826924 : if( previousVerticesMemento )
530 : 826814 : pd.recreate_vertices_memento( previousVerticesMemento, err );
531 : : else
532 [ - + ][ # # ]: 826924 : previousVerticesMemento = pd.create_vertices_memento( err );MSQ_ERRRTN( err );
[ - + ]
533 : : }
534 : :
535 [ + + ]: 1659896 : if( totalFlag & UNTANGLED_MESH )
536 : : {
537 : 57185 : patchInvertedCount = count_inverted( pd, err );
538 : : // MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << " o Num Patch Inverted: " << " " <<
539 : : // patchInvertedCount << std::endl;MSQ_ERRRTN( err );
540 : : }
541 : : }
542 : :
543 : 832 : void TerminationCriterion::accumulate_inner( PatchData& pd, OFEvaluator& of_eval, MsqError& err )
544 : : {
545 : 832 : double of_value = 0;
546 : :
547 [ - + ]: 832 : if( terminationCriterionFlag & GRAD_FLAGS )
548 : : {
549 [ # # ][ # # ]: 0 : mGrad.resize( pd.num_free_vertices() );
550 [ # # ][ # # ]: 832 : bool b = of_eval.evaluate( pd, of_value, mGrad, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
551 [ # # ]: 0 : if( !b )
552 : : {
553 : : MSQ_SETERR( err )
554 [ # # ][ # # ]: 0 : ( "Initial patch is invalid for gradient compuation.", MsqError::INVALID_MESH );
555 : 0 : return;
556 : : }
557 : : }
558 [ + + ]: 832 : else if( terminationCriterionFlag & OF_FLAGS )
559 : : {
560 [ + - ][ + - ]: 4 : bool b = of_eval.evaluate( pd, of_value, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
561 [ - + ]: 4 : if( !b )
562 : : {
563 : : MSQ_SETERR( err )
564 [ # # ][ # # ]: 0 : ( "Invalid patch passed to TerminationCriterion.", MsqError::INVALID_MESH );
565 : 4 : return;
566 : : }
567 : : }
568 : :
569 [ - + ][ + - ]: 832 : accumulate_inner( pd, of_value, mGrad.empty() ? 0 : arrptr( mGrad ), err );MSQ_CHKERR( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ]
570 : : }
571 : :
572 : 60670 : void TerminationCriterion::accumulate_inner( PatchData& pd, double of_value, Vector3D* grad_array, MsqError& err )
573 : : {
574 : : // if terminating on the norm of the gradient
575 : : // currentGradL2NormSquared = HUGE_VAL;
576 [ + + ]: 60670 : if( terminationCriterionFlag & ( GRADIENT_L2_NORM_ABSOLUTE | GRADIENT_L2_NORM_RELATIVE ) )
577 : : {
578 : 13 : currentGradL2NormSquared = length_squared( grad_array, pd.num_free_vertices() ); // get the L2 norm
579 : 13 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o Info -- gradient L2 norm: "
580 : : << " " << RPM( std::sqrt( currentGradL2NormSquared ) ) << std::endl;
581 : : }
582 : : // currentGradInfNorm = 10e6;
583 [ - + ]: 60670 : if( terminationCriterionFlag & ( GRADIENT_INF_NORM_ABSOLUTE | GRADIENT_INF_NORM_RELATIVE ) )
584 : : {
585 : 0 : currentGradInfNorm = Linf( grad_array, pd.num_free_vertices() ); // get the Linf norm
586 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o Info -- gradient Inf norm: "
587 : : << " " << RPM( currentGradInfNorm ) << std::endl;
588 : : }
589 : :
590 [ - + ]: 60670 : if( terminationCriterionFlag & VERTEX_MOVEMENT_RELATIVE )
591 : : {
592 [ # # ][ # # ]: 0 : maxSquaredInitialMovement = pd.get_max_vertex_movement_squared( initialVerticesMemento, err );MSQ_ERRRTN( err );
[ # # ]
593 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o Info -- max initial vertex movement: "
594 : : << " " << RPM( maxSquaredInitialMovement ) << std::endl;
595 : : }
596 : :
597 : 60670 : previousOFValue = currentOFValue;
598 : 60670 : currentOFValue = of_value;
599 [ + + ]: 60670 : if( terminationCriterionFlag & OF_FLAGS )
600 : : {
601 : 26 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o Info -- OF Value: "
602 : : << " " << RPM( of_value ) << " iterationCounter= " << iterationCounter
603 : : << std::endl;
604 : : }
605 [ + + ]: 60644 : else if( grad_array )
606 : : {
607 : 40784 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o OF Value: "
608 : : << " " << RPM( of_value ) << " iterationCounter= "
609 : : << iterationCounter
610 : : //<< " terminationCriterionFlag= " <<
611 : : // terminationCriterionFlag << " OF_FLAGS = " << OF_FLAGS
612 : : << std::endl;
613 : : }
614 : :
615 : 60670 : ++iterationCounter;
616 [ + + ]: 60670 : if( timeStepFileType ) write_timestep( pd, grad_array, err );
617 : :
618 [ - + ]: 60670 : if( plotFile.is_open() )
619 : 0 : plotFile << iterationCounter << '\t' << mTimer.since_birth() << '\t' << of_value << '\t'
620 : 0 : << std::sqrt( currentGradL2NormSquared ) << '\t' << currentGradInfNorm << '\t'
621 [ # # ]: 0 : << ( maxSquaredMovement > 0.0 ? std::sqrt( maxSquaredMovement ) : 0.0 ) << '\t' << globalInvertedCount
622 : 0 : << std::endl;
623 : : }
624 : :
625 : 832 : void TerminationCriterion::accumulate_outer( Mesh* mesh, MeshDomain* domain, OFEvaluator& of_eval,
626 : : const Settings* settings, MsqError& err )
627 : : {
628 [ + - ]: 832 : PatchData global_patch;
629 [ + - ][ + - ]: 832 : if( settings ) global_patch.attach_settings( settings );
630 : :
631 : : // if we need to fill out the global patch data object.
632 [ + + ][ - + ]: 832 : if( ( terminationCriterionFlag & ( GRAD_FLAGS | OF_FLAGS | VERTEX_MOVEMENT_RELATIVE ) ) || timeStepFileType )
633 : : {
634 [ + - ]: 4 : global_patch.set_mesh( mesh );
635 [ + - ]: 4 : global_patch.set_domain( domain );
636 [ + - ][ + - ]: 4 : global_patch.fill_global_patch( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
637 : : }
638 : :
639 [ + - ][ + - ]: 832 : accumulate_inner( global_patch, of_eval, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
[ + - ]
640 : : }
641 : :
642 : 869924 : void TerminationCriterion::accumulate_patch( PatchData& pd, MsqError& err )
643 : : {
644 [ + + ]: 869924 : if( terminationCriterionFlag & MOVEMENT_FLAGS )
645 : : {
646 : 980 : double patch_max_dist = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
647 [ + - ]: 980 : if( patch_max_dist > maxSquaredMovement ) maxSquaredMovement = patch_max_dist;
648 [ - + ][ # # ]: 980 : pd.recreate_vertices_memento( previousVerticesMemento, err );MSQ_ERRRTN( err );
[ - + ]
649 : : }
650 : :
651 : : // if terminating on bounded vertex movement (a bounding box for the mesh)
652 [ - + ]: 869924 : if( terminationCriterionFlag & BOUNDED_VERTEX_MOVEMENT )
653 : : {
654 : 0 : const MsqVertex* vert = pd.get_vertex_array( err );
655 : 0 : int num_vert = pd.num_free_vertices();
656 : 0 : int i = 0;
657 : : // for each vertex
658 [ # # ]: 0 : for( i = 0; i < num_vert; ++i )
659 : : {
660 : : // if any of the coordinates are greater than eps
661 [ # # ]: 0 : if( ( vert[i][0] > boundedVertexMovementEps ) || ( vert[i][1] > boundedVertexMovementEps ) ||
[ # # # # ]
[ # # ]
662 : 0 : ( vert[i][2] > boundedVertexMovementEps ) )
663 : 0 : { ++vertexMovementExceedsBound; }
664 : : }
665 : : }
666 : :
667 [ + + ]: 869924 : if( ( terminationCriterionFlag | cullingMethodFlag ) & UNTANGLED_MESH )
668 : : {
669 : 57185 : size_t new_count = count_inverted( pd, err );
670 : : // be careful here because size_t is unsigned
671 : 57185 : globalInvertedCount += new_count;
672 : 57185 : globalInvertedCount -= patchInvertedCount;
673 : 57185 : patchInvertedCount = new_count;
674 : : // if (innerOuterType==TYPE_OUTER)
675 : : // MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << " o Num Patch Inverted: " << " " <<
676 : : // patchInvertedCount << " globalInvertedCount= " << globalInvertedCount << std::endl;
677 : :
678 [ - + ][ # # ]: 57185 : MSQ_ERRRTN( err );
[ - + ]
679 : : }
680 : : }
681 : :
682 : : /*! This function evaluates the needed information and then evaluates
683 : : the termination criteria. If any of the selected criteria are satisfied,
684 : : the function returns true. Otherwise, the function returns false.
685 : : */
686 : 939474 : bool TerminationCriterion::terminate()
687 : : {
688 : 939474 : bool return_flag = false;
689 : : // std::cout<<"\nInside terminate(pd,of,err): flag = "<<terminationCriterionFlag << std::endl;
690 : 939474 : int type = 0;
691 : :
692 : : // First check for an interrupt signal
693 [ - + ]: 939474 : if( MsqInterrupt::interrupt() )
694 : : {
695 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- INTERRUPTED"
696 : : << " " << std::endl;
697 : 0 : return true;
698 : : }
699 : :
700 : : // if terminating on numbering of inner iterations
701 [ + + ][ + + ]: 939474 : if( NUMBER_OF_ITERATES & terminationCriterionFlag && iterationCounter >= iterationBound )
702 : : {
703 : 39616 : return_flag = true;
704 : 39616 : type = 1;
705 : 39616 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- Reached "
706 : : << " " << iterationBound << " iterations." << std::endl;
707 : : }
708 : :
709 [ + + ][ - + ]: 939474 : if( CPU_TIME & terminationCriterionFlag && mTimer.since_birth() >= timeBound )
[ - + ]
710 : : {
711 : 0 : return_flag = true;
712 : 0 : type = 2;
713 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- Exceeded CPU time. "
714 : : << " " << std::endl;
715 : : }
716 : :
717 [ + + ][ + + ]: 939474 : if( MOVEMENT_FLAGS & terminationCriterionFlag && maxSquaredMovement >= 0.0 )
718 : : {
719 : 980 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o Info -- Maximum vertex movement: "
720 : : << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
721 : :
722 [ + + ][ + + ]: 980 : if( VERTEX_MOVEMENT_ABSOLUTE & terminationCriterionFlag && maxSquaredMovement <= vertexMovementAbsoluteEps )
723 : : {
724 : 79 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- VERTEX_MOVEMENT_ABSOLUTE: "
725 : : << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
726 : 79 : return_flag = true;
727 : 79 : type = 3;
728 : : }
729 : :
730 [ - + ][ # # ]: 980 : if( VERTEX_MOVEMENT_RELATIVE & terminationCriterionFlag &&
731 : 0 : maxSquaredMovement <= vertexMovementRelativeEps * maxSquaredInitialMovement )
732 : : {
733 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- VERTEX_MOVEMENT_RELATIVE: "
734 : : << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
735 : 0 : return_flag = true;
736 : 0 : type = 4;
737 : : }
738 : :
739 [ + + ]: 980 : if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & terminationCriterionFlag )
740 : : {
741 [ - + ]: 58 : assert( vertexMovementAbsoluteAvgEdge > -1e-12 ); // make sure value actually got calculated
742 [ + + ]: 58 : if( maxSquaredMovement <= vertexMovementAbsoluteAvgEdge )
743 : : {
744 : 3 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- VERTEX_MOVEMENT_ABS_EDGE_LENGTH: "
745 : : << " " << RPM( sqrt( maxSquaredMovement ) ) << std::endl;
746 : 3 : return_flag = true;
747 : 3 : type = 5;
748 : : }
749 : : }
750 : : }
751 : :
752 [ + + ][ + + ]: 939474 : if( GRADIENT_L2_NORM_ABSOLUTE & terminationCriterionFlag &&
753 : 18 : currentGradL2NormSquared <= gradL2NormAbsoluteEpsSquared )
754 : : {
755 : 1 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- GRADIENT_L2_NORM_ABSOLUTE: "
756 : : << " " << RPM( currentGradL2NormSquared ) << std::endl;
757 : 1 : return_flag = true;
758 : 1 : type = 6;
759 : : }
760 : :
761 [ - + ][ # # ]: 939474 : if( GRADIENT_INF_NORM_ABSOLUTE & terminationCriterionFlag && currentGradInfNorm <= gradInfNormAbsoluteEps )
762 : : {
763 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- GRADIENT_INF_NORM_ABSOLUTE: "
764 : : << " " << RPM( currentGradInfNorm ) << std::endl;
765 : 0 : return_flag = true;
766 : 0 : type = 7;
767 : : }
768 : :
769 [ - + ][ # # ]: 939474 : if( GRADIENT_L2_NORM_RELATIVE & terminationCriterionFlag &&
770 : 0 : currentGradL2NormSquared <= ( gradL2NormRelativeEpsSquared * initialGradL2NormSquared ) )
771 : : {
772 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- GRADIENT_L2_NORM_RELATIVE: "
773 : : << " " << RPM( currentGradL2NormSquared ) << std::endl;
774 : 0 : return_flag = true;
775 : 0 : type = 8;
776 : : }
777 : :
778 [ - + ][ # # ]: 939474 : if( GRADIENT_INF_NORM_RELATIVE & terminationCriterionFlag &&
779 : 0 : currentGradInfNorm <= ( gradInfNormRelativeEps * initialGradInfNorm ) )
780 : : {
781 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- GRADIENT_INF_NORM_RELATIVE: "
782 : : << " " << RPM( currentGradInfNorm ) << std::endl;
783 : 0 : return_flag = true;
784 : 0 : type = 9;
785 : : }
786 : : // Quality Improvement and Successive Improvements are below.
787 : : // The relative forms are only valid after the first iteration.
788 [ + + ][ + + ]: 939474 : if( ( QUALITY_IMPROVEMENT_ABSOLUTE & terminationCriterionFlag ) && currentOFValue <= qualityImprovementAbsoluteEps )
789 : : {
790 : 3 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- QUALITY_IMPROVEMENT_ABSOLUTE: "
791 : : << " " << RPM( currentOFValue ) << std::endl;
792 : 3 : return_flag = true;
793 : 3 : type = 10;
794 : : }
795 : :
796 : : // only valid if an iteration has occurred, see above.
797 [ + + ]: 939474 : if( iterationCounter > 0 )
798 : : {
799 [ - + ][ # # ]: 51211 : if( SUCCESSIVE_IMPROVEMENTS_ABSOLUTE & terminationCriterionFlag &&
800 : 0 : ( previousOFValue - currentOFValue ) <= successiveImprovementsAbsoluteEps )
801 : : {
802 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel )
803 : : << par_string() << " o TermCrit -- SUCCESSIVE_IMPROVEMENTS_ABSOLUTE: previousOFValue= "
804 : : << " " << previousOFValue << " currentOFValue= " << RPM( currentOFValue )
805 : : << " successiveImprovementsAbsoluteEps= " << successiveImprovementsAbsoluteEps << std::endl;
806 : 0 : return_flag = true;
807 : 0 : type = 11;
808 : : }
809 [ + + ][ + + ]: 51211 : if( QUALITY_IMPROVEMENT_RELATIVE & terminationCriterionFlag &&
810 : 15 : ( currentOFValue - lowerOFBound ) <= qualityImprovementRelativeEps * ( initialOFValue - lowerOFBound ) )
811 : : {
812 : 2 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- QUALITY_IMPROVEMENT_RELATIVE: "
813 : : << " " << std::endl;
814 : 2 : return_flag = true;
815 : 2 : type = 12;
816 : : }
817 [ + + ][ + + ]: 51211 : if( SUCCESSIVE_IMPROVEMENTS_RELATIVE & terminationCriterionFlag &&
818 : 10 : ( previousOFValue - currentOFValue ) <=
819 : 10 : successiveImprovementsRelativeEps * ( initialOFValue - currentOFValue ) )
820 : : {
821 : 1 : MSQ_DBGOUT_P0_ONLY( debugLevel )
822 : : << par_string() << " o TermCrit -- SUCCESSIVE_IMPROVEMENTS_RELATIVE: previousOFValue= "
823 : : << " " << previousOFValue << " currentOFValue= " << RPM( currentOFValue )
824 : : << " successiveImprovementsRelativeEps= " << successiveImprovementsRelativeEps << std::endl;
825 : 1 : return_flag = true;
826 : 1 : type = 13;
827 : : }
828 : : }
829 : :
830 [ - + ][ # # ]: 939474 : if( BOUNDED_VERTEX_MOVEMENT & terminationCriterionFlag && vertexMovementExceedsBound )
831 : : {
832 : 0 : return_flag = true;
833 : 0 : type = 14;
834 : 0 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- "
835 : : << " " << vertexMovementExceedsBound << " vertices out of bounds."
836 : : << std::endl;
837 : : }
838 : :
839 [ + + ]: 939474 : if( UNTANGLED_MESH & terminationCriterionFlag )
840 : : {
841 : 344 : if( innerOuterType == TYPE_OUTER )
842 : : {
843 : : // MSQ_DBGOUT_P0_ONLY(debugLevel)
844 : : MSQ_DBGOUT( debugLevel ) << par_string() << " o Num Inverted: "
845 : : << " " << globalInvertedCount << std::endl;
846 : : }
847 [ + + ]: 344 : if( !globalInvertedCount )
848 : : {
849 : 18 : MSQ_DBGOUT_P0_ONLY( debugLevel ) << par_string() << " o TermCrit -- UNTANGLED_MESH: "
850 : : << " " << std::endl;
851 : 18 : return_flag = true;
852 : 18 : type = 15;
853 : : }
854 : : }
855 : :
856 : : // clear this value at the end of each iteration
857 : 939474 : vertexMovementExceedsBound = 0;
858 : 939474 : maxSquaredMovement = -1.0;
859 : :
860 [ + + ][ + + ]: 939474 : if( timeStepFileType == GNUPLOT && return_flag )
861 : : {
862 [ + - ]: 1 : MsqError err;
863 [ + - ]: 939474 : MeshWriter::write_gnuplot_overlay( iterationCounter, timeStepFileName.c_str(), err );
864 : : }
865 : :
866 : : if( 0 && return_flag && MSQ_DBG( 2 ) )
867 : : std::cout << "P[" << get_parallel_rank() << "] tmp TerminationCriterion::terminate: " << moniker
868 : : << " return_flag= " << return_flag << " type= " << type
869 : : << " terminationCriterionFlag= " << terminationCriterionFlag << " debugLevel= " << debugLevel
870 : : << std::endl;
871 : :
872 : : // if none of the criteria were satisfied
873 : 939474 : return return_flag;
874 : : }
875 : :
876 : 256 : bool TerminationCriterion::criterion_is_set()
877 : : {
878 [ + + ]: 256 : if( !terminationCriterionFlag )
879 : 11 : return false;
880 : : else
881 : 245 : return true;
882 : : }
883 : :
884 : : /*!This function checks the culling method criterion supplied to the object
885 : : by the user. If the user does not supply a culling method criterion,
886 : : the default criterion is NONE, and in that case, no culling is performed.
887 : : If the culling method criterion is satisfied, the interior vertices
888 : : of the given patch are flagged as soft_fixed. Otherwise, the soft_fixed
889 : : flag is removed from each of the vertices in the patch (interior and
890 : : boundary vertices). Also, if the criterion was satisfied, then the
891 : : function returns true. Otherwise, the function returns false.
892 : : */
893 : 829946 : bool TerminationCriterion::cull_vertices( PatchData& pd, OFEvaluator& of_eval, MsqError& err )
894 : : {
895 : : // PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
896 : :
897 : : // cull_bool will be changed to true if the criterion is satisfied
898 : 829946 : bool b, cull_bool = false;
899 : : double prev_m, init_m;
900 [ + - - + : 829946 : switch( cullingMethodFlag )
- - - ]
901 : : {
902 : : // if no culling is requested, always return false
903 : : case NONE:
904 : 3833 : return cull_bool;
905 : : // if culling on quality improvement absolute
906 : : case QUALITY_IMPROVEMENT_ABSOLUTE:
907 : : // get objective function value
908 : 0 : b = of_eval.evaluate( pd, currentOFValue, err );
909 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) ) return false;
[ # # ]
910 [ # # ]: 0 : if( !b )
911 : : {
912 [ # # ]: 0 : MSQ_SETERR( err )( MsqError::INVALID_MESH );
913 : 0 : return false;
914 : : }
915 : : // if the improvement was enough, cull
916 [ # # ]: 0 : if( currentOFValue <= cullingEps ) { cull_bool = true; }
917 : : // PRINT_INFO("\ncurrentOFValue = %f, bool = %i\n",currentOFValue,cull_bool);
918 : :
919 : 0 : break;
920 : : // if culing on quality improvement relative
921 : : case QUALITY_IMPROVEMENT_RELATIVE:
922 : : // get objective function value
923 : 0 : b = of_eval.evaluate( pd, currentOFValue, err );
924 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) ) return false;
[ # # ]
925 [ # # ]: 0 : if( !b )
926 : : {
927 [ # # ]: 0 : MSQ_SETERR( err )( MsqError::INVALID_MESH );
928 : 0 : return false;
929 : : }
930 : : // if the improvement was enough, cull
931 [ # # ]: 0 : if( ( currentOFValue - lowerOFBound ) <= ( cullingEps * ( initialOFValue - lowerOFBound ) ) )
932 : 0 : { cull_bool = true; }
933 : 0 : break;
934 : : // if culling on vertex movement absolute
935 : : case VERTEX_MOVEMENT_ABSOLUTE:
936 : : case VERTEX_MOVEMENT_ABS_EDGE_LENGTH:
937 : : // if movement was enough, cull
938 : 826113 : prev_m = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
939 [ - + ][ # # ]: 826113 : MSQ_ERRZERO( err );
[ - + ]
940 [ + + ]: 826113 : if( prev_m <= cullingEps * cullingEps ) { cull_bool = true; }
941 : :
942 : 826113 : break;
943 : : // if culling on vertex movement relative
944 : : case VERTEX_MOVEMENT_RELATIVE:
945 : : // if movement was small enough, cull
946 : 0 : prev_m = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
947 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ]
948 : 0 : init_m = pd.get_max_vertex_movement_squared( initialVerticesMemento, err );
949 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ]
950 [ # # ]: 0 : if( prev_m <= ( cullingEps * cullingEps * init_m ) ) { cull_bool = true; }
951 : 0 : break;
952 : : case UNTANGLED_MESH:
953 [ # # ]: 0 : if( !patchInvertedCount ) cull_bool = true;
954 : 0 : break;
955 : : default:
956 : : MSQ_SETERR( err )
957 [ # # ]: 0 : ( "Requested culling method not yet implemented.", MsqError::NOT_IMPLEMENTED );
958 : 0 : return false;
959 : : };
960 : : // Now actually have patch data cull vertices
961 [ + + ]: 826113 : if( cull_bool )
962 : : {
963 : 49295 : pd.set_free_vertices_soft_fixed( err );
964 [ - + ][ # # ]: 49295 : MSQ_ERRZERO( err );
[ - + ]
965 : : }
966 : : else
967 : : {
968 : 776818 : pd.set_all_vertices_soft_free( err );
969 [ - + ][ # # ]: 776818 : MSQ_ERRZERO( err );
[ - + ]
970 : : }
971 : 829946 : return cull_bool;
972 : : }
973 : :
974 : : /*!This function is activated when cullingGlobalPatch is true. It supplies
975 : : cull_vertices with a single vertex-based patch at a time. If the patch
976 : : satisfies the culling criterion, it's free vertices are then soft-fixed.
977 : : */
978 : 0 : bool TerminationCriterion::cull_vertices_global( PatchData& global_patch, Mesh* mesh, MeshDomain* domain,
979 : : const Settings* settings, OFEvaluator& of_eval, MsqError& err )
980 : : {
981 [ # # ]: 0 : if( !cullingGlobalPatch ) return false;
982 : :
983 : : // PRINT_INFO("CULLING_METHOD FLAG = %i",cullingMethodFlag);
984 : :
985 : : // cull_bool will be changed to true if the criterion is satisfied
986 : 0 : bool cull_bool = false;
987 : :
988 [ # # ]: 0 : std::vector< Mesh::VertexHandle > mesh_vertices;
989 : : // std::vector<Mesh::VertexHandle> patch_vertices;
990 : : // std::vector<Mesh::ElementHandle> patch_elements;
991 : : // std::vector<Mesh::VertexHandle> fixed_vertices;
992 : : // std::vector<Mesh::VertexHandle> free_vertices;
993 : :
994 : : // FIXME, verify global_patch is a global patch... how, is this right?
995 [ # # ]: 0 : mesh->get_all_vertices( mesh_vertices, err );
996 : 0 : size_t mesh_num_nodes = mesh_vertices.size();
997 [ # # ]: 0 : size_t global_patch_num_nodes = global_patch.num_nodes();
998 : : if( 0 )
999 : : std::cout << "tmp srk mesh_num_nodes= " << mesh_num_nodes
1000 : : << " global_patch_num_nodes= " << global_patch_num_nodes << std::endl;
1001 [ # # ]: 0 : if( mesh_num_nodes != global_patch_num_nodes )
1002 : : {
1003 [ # # ][ # # ]: 0 : std::cout << "tmp srk cull_vertices_global found non global patch" << std::endl;
1004 : 0 : exit( 123 );
1005 : : return false;
1006 : : }
1007 [ # # ]: 0 : PatchData patch;
1008 [ # # ]: 0 : patch.set_mesh( (Mesh*)mesh );
1009 [ # # ]: 0 : patch.set_domain( domain );
1010 [ # # ]: 0 : patch.attach_settings( settings );
1011 : :
1012 : : // const MsqVertex* global_patch_vertex_array = global_patch.get_vertex_array( err );
1013 [ # # ]: 0 : Mesh::VertexHandle* global_patch_vertex_handles = global_patch.get_vertex_handles_array();
1014 : :
1015 : 0 : int num_culled = 0;
1016 [ # # ]: 0 : for( unsigned iv = 0; iv < global_patch_num_nodes; iv++ )
1017 : : {
1018 : : // form a patch for this vertex; if it is culled, set it to be soft fixed
1019 : 0 : Mesh::VertexHandle vert = global_patch_vertex_handles[iv];
1020 [ # # ]: 0 : std::vector< Mesh::ElementHandle > elements;
1021 [ # # ]: 0 : std::vector< size_t > offsets;
1022 [ # # ]: 0 : mesh->vertices_get_attached_elements( &vert, 1, elements, offsets, err );
1023 : :
1024 [ # # ]: 0 : std::set< Mesh::VertexHandle > patch_free_vertices_set;
1025 : :
1026 [ # # ]: 0 : for( unsigned ie = 0; ie < elements.size(); ie++ )
1027 : : {
1028 [ # # ]: 0 : std::vector< Mesh::VertexHandle > vert_handles;
1029 [ # # ]: 0 : std::vector< size_t > v_offsets;
1030 [ # # ][ # # ]: 0 : mesh->elements_get_attached_vertices( &elements[ie], 1, vert_handles, v_offsets, err );
1031 [ # # ]: 0 : for( unsigned jv = 0; jv < vert_handles.size(); jv++ )
1032 : : {
1033 : : unsigned char bt;
1034 [ # # ][ # # ]: 0 : mesh->vertex_get_byte( vert_handles[jv], &bt, err );
1035 [ # # ]: 0 : MsqVertex v;
1036 [ # # ]: 0 : v.set_flags( bt );
1037 [ # # ][ # # ]: 0 : if( v.is_free_vertex() ) patch_free_vertices_set.insert( vert_handles[jv] );
[ # # ][ # # ]
1038 : : }
1039 : 0 : }
1040 : :
1041 : : std::vector< Mesh::VertexHandle > patch_free_vertices_vector( patch_free_vertices_set.begin(),
1042 [ # # ]: 0 : patch_free_vertices_set.end() );
1043 : : // std::vector<unsigned char> byte_vector(patch_vertices_vector.size());
1044 : : // mesh->vertices_get_byte(&vert_handles[0], &byte_vector[0], vert_handles.size(), err);
1045 : :
1046 [ # # ]: 0 : patch.set_mesh_entities( elements, patch_free_vertices_vector, err );
1047 [ # # ][ # # ]: 0 : if( cull_vertices( patch, of_eval, err ) )
1048 : : {
1049 : : // std::cout << "tmp srk cull_vertices_global found culled patch" << std::endl;
1050 [ # # ]: 0 : Mesh::VertexHandle* patch_vertex_handles = patch.get_vertex_handles_array();
1051 [ # # ]: 0 : const MsqVertex* patch_vertex_array = patch.get_vertex_array( err );
1052 [ # # ][ # # ]: 0 : for( unsigned jv = 0; jv < patch.num_nodes(); jv++ )
1053 : : {
1054 [ # # ]: 0 : if( patch_vertex_handles[jv] == global_patch_vertex_handles[iv] )
1055 : : {
1056 [ # # ][ # # ]: 0 : if( patch_vertex_array[jv].is_flag_set( MsqVertex::MSQ_CULLED ) )
1057 : : {
1058 [ # # ]: 0 : global_patch.set_vertex_culled( iv );
1059 : 0 : ++num_culled;
1060 : 0 : cull_bool = true;
1061 : : // std::cout << "tmp srk cull_vertices_global found culled vertex" <<
1062 : : // std::endl;
1063 : : }
1064 : : }
1065 : : }
1066 : : }
1067 : 0 : }
1068 : : if( 0 )
1069 : : std::cout << "tmp srk cull_vertices_global found " << num_culled << " culled vertices out of "
1070 : : << global_patch_num_nodes << std::endl;
1071 : :
1072 : 0 : return cull_bool;
1073 : : }
1074 : :
1075 : 114390 : size_t TerminationCriterion::count_inverted( PatchData& pd, MsqError& err )
1076 : : {
1077 [ + - ]: 114390 : size_t num_elem = pd.num_elements();
1078 : 114390 : size_t count = 0;
1079 : : int inverted, samples;
1080 [ + + ]: 582630 : for( size_t i = 0; i < num_elem; i++ )
1081 : : {
1082 [ + - ][ + - ]: 468240 : pd.element_by_index( i ).check_element_orientation( pd, inverted, samples, err );
1083 [ + + ]: 468240 : if( inverted ) ++count;
1084 : : }
1085 : 114390 : return count;
1086 : : }
1087 : :
1088 : : /*!
1089 : : Currently this only deletes the memento of the vertex positions and the
1090 : : mGrad vector if neccessary.
1091 : : When culling, we remove the soft fixed flags from all of the vertices.
1092 : : */
1093 : 256 : void TerminationCriterion::cleanup( Mesh*, MeshDomain*, MsqError& )
1094 : : {
1095 [ + + ]: 256 : delete previousVerticesMemento;
1096 [ - + ]: 256 : delete initialVerticesMemento;
1097 : 256 : previousVerticesMemento = 0;
1098 : 256 : initialVerticesMemento = 0;
1099 : 256 : }
1100 : :
1101 : 11 : void TerminationCriterion::write_timestep( PatchData& pd, const Vector3D* gradient, MsqError& err )
1102 : : {
1103 [ + - ]: 11 : std::ostringstream str;
1104 [ - + ]: 11 : if( timeStepFileType == VTK )
1105 : : {
1106 [ # # ][ # # ]: 0 : str << timeStepFileName << '_' << iterationCounter << ".vtk";
[ # # ][ # # ]
1107 [ # # ][ # # ]: 0 : MeshWriter::write_vtk( pd, str.str().c_str(), err, gradient );
1108 : : }
1109 [ + - ]: 11 : else if( timeStepFileType == GNUPLOT )
1110 : : {
1111 [ + - ][ + - ]: 11 : str << timeStepFileName << '.' << iterationCounter;
[ + - ]
1112 [ + - ][ + - ]: 11 : MeshWriter::write_gnuplot( pd.get_mesh(), str.str().c_str(), err );
[ + - ]
1113 : 11 : }
1114 : 11 : }
1115 : :
1116 : 0 : void TerminationCriterion::write_iterations( const char* filename, MsqError& err )
1117 : : {
1118 [ # # ]: 0 : if( filename )
1119 : : {
1120 : 0 : plotFile.open( filename, std::ios::trunc );
1121 [ # # ]: 0 : if( !plotFile ) MSQ_SETERR( err )
1122 [ # # ]: 0 : ( MsqError::FILE_ACCESS, "Failed to open plot data file: '%s'", filename );
1123 : : }
1124 : : else
1125 : : {
1126 : 0 : plotFile.close();
1127 : : }
1128 : 0 : }
1129 : :
1130 : 260 : void TerminationCriterion::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings*, MsqError& err )
1131 : : {
1132 [ + + ]: 260 : if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & ( terminationCriterionFlag | cullingMethodFlag ) )
1133 : : {
1134 [ + - ]: 10 : Mesh* mesh = mesh_and_domain->get_mesh();
1135 [ + - ]: 10 : MeshUtil tool( mesh );
1136 [ + - ]: 10 : SimpleStats stats;
1137 [ + - ][ + - ]: 270 : tool.edge_length_distribution( stats, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
1138 [ + - ][ + - ]: 10 : double limit = vertexMovementAvgBeta * ( stats.average() - stats.standard_deviation() );
1139 : : // we actually calculate the square of the length
1140 : 10 : vertexMovementAbsoluteAvgEdge = limit * limit;
1141 [ + + ][ + - ]: 10 : if( VERTEX_MOVEMENT_ABS_EDGE_LENGTH & cullingMethodFlag ) cullingEps = limit;
1142 : : }
1143 : : }
1144 : :
1145 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
