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