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 : : /*!
28 : : \file VertexMover.cpp
29 : : \brief
30 : :
31 : : The VertexMover Class is the base class for all the smoothing algorythms
32 : :
33 : : \author Thomas Leurent
34 : : \date 2002-01-17
35 : : */
36 : :
37 : : #include "VertexMover.hpp"
38 : : #include "NonGradient.hpp"
39 : : #include "QualityMetric.hpp"
40 : : #include "TQualityMetric.hpp"
41 : : #include "TMetricBarrier.hpp"
42 : : #include "AWMetricBarrier.hpp"
43 : : #include "ElementMaxQM.hpp"
44 : : #include "ElementAvgQM.hpp"
45 : : #include "ElementPMeanP.hpp"
46 : : #include "PlanarDomain.hpp"
47 : : #include "ObjectiveFunctionTemplate.hpp"
48 : : #include "MaxTemplate.hpp"
49 : : #include "MsqTimer.hpp"
50 : : #include "MsqDebug.hpp"
51 : : #include "PatchSet.hpp"
52 : : #include "PatchData.hpp"
53 : : #include "ParallelHelperInterface.hpp"
54 : : #include <algorithm>
55 : :
56 : : namespace MBMesquite
57 : : {
58 : :
59 : : extern int get_parallel_rank();
60 : : extern int get_parallel_size();
61 : : extern void parallel_barrier();
62 : :
63 [ + - ]: 163 : VertexMover::VertexMover( ObjectiveFunction* OF ) : QualityImprover(), objFuncEval( OF ), jacobiOpt( false ) {}
64 : :
65 [ - + ]: 312 : VertexMover::~VertexMover() {}
66 : :
67 : : /*
68 : :
69 : : +-----------+
70 : : |Reset Outer|
71 : : |Criterion |
72 : : +-----------+
73 : : |
74 : : V
75 : : +
76 : : / \
77 : : /Outer\ YES
78 : : +--> <Criterion>-----> DONE
79 : : | \Done?/
80 : : | \ /
81 : : | +
82 : : | |NO
83 : : | V
84 : : | +-----------+
85 : : 1 |Reset Mesh |
86 : : | | Iteration |
87 : : | +-----------+
88 : : | |
89 : : | V
90 : : | +
91 : : | / \
92 : : | NO /Next\
93 : : +-----<Patch > <-----+
94 : : \ / |
95 : : \ / |
96 : : + |
97 : : YES| |
98 : : V |
99 : : +-----------+ |
100 : : |Reset Inner| |
101 : : |Criterion | 2
102 : : +-----------+ |
103 : : | |
104 : : V |
105 : : + |
106 : : / \ |
107 : : /Inner\ YES |
108 : : +--> <Criterion>-----+ --------------+
109 : : | \Done?/ |
110 : : | \ / |
111 : : | + |
112 : : | |NO |
113 : : | V Inside
114 : : 3 +-----------+ Smoother
115 : : | | Smooth | |
116 : : | | Patch | |
117 : : | +-----------+ |
118 : : | | |
119 : : ----------+ --------------+
120 : :
121 : : */
122 : :
123 : : /*! \brief Improves the quality of the MeshSet, calling some
124 : : methods specified in a class derived from VertexMover
125 : :
126 : : \param const MeshSet &: this MeshSet is looped over. Only the
127 : : mutable data members are changed (such as currentVertexInd).
128 : : */
129 : 130 : double VertexMover::loop_over_mesh( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err )
130 : : {
131 [ + - ]: 130 : Mesh* mesh = mesh_and_domain->get_mesh();
132 [ + - ]: 130 : MeshDomain* domain = mesh_and_domain->get_domain();
133 : :
134 : 130 : TagHandle coord_tag = 0; // store uncommitted coords for jacobi optimization
135 : 130 : TagHandle* coord_tag_ptr = 0;
136 : :
137 : 130 : TerminationCriterion* outer_crit = 0;
138 : 130 : TerminationCriterion* inner_crit = 0;
139 : :
140 : : // Clear culling flag, set hard fixed flag, etc on all vertices
141 [ + - ]: 130 : initialize_vertex_byte( mesh_and_domain, settings, err );
142 [ + - ][ - + ]: 130 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
143 : :
144 : : // Get the patch data to use for the first iteration
145 [ + - ]: 130 : OFEvaluator& obj_func = get_objective_function_evaluator();
146 : :
147 : : // Check for impromer use of MaxTemplate
148 [ + - ][ + + ]: 130 : ObjectiveFunctionTemplate* maxt_ptr = dynamic_cast< MaxTemplate* >( obj_func.get_objective_function() );
149 [ + + ]: 130 : if( maxt_ptr )
150 : : {
151 [ - + ]: 2 : QualityImprover* ngqi_ptr = dynamic_cast< NonGradient* >( this );
152 [ - + ]: 2 : if( !ngqi_ptr )
153 [ # # ][ # # ]: 0 : std::cout << "Warning: MaxTemplate results in non-differentiable objective function." << std::endl
154 [ # # ][ # # ]: 0 : << " Therefore, it is best to use the NonGradient solver. Other Mesquite" << std::endl
155 [ # # ][ # # ]: 2 : << " solvers require derivative information." << std::endl;
156 : : }
157 : :
158 [ + - ]: 130 : PatchData patch;
159 [ + - ]: 130 : patch.set_mesh( mesh );
160 [ + - ]: 130 : patch.set_domain( domain );
161 [ + - ][ + - ]: 130 : if( settings ) patch.attach_settings( settings );
162 : 130 : bool one_patch = false, inner_crit_terminated, all_culled;
163 [ + - ]: 260 : std::vector< Mesh::VertexHandle > patch_vertices;
164 [ + - ]: 260 : std::vector< Mesh::ElementHandle > patch_elements;
165 : : bool valid;
166 : :
167 [ + - ]: 130 : PatchSet* patch_set = get_patch_set();
168 [ - + ]: 130 : if( !patch_set )
169 : : {
170 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No PatchSet for QualityImprover!", MsqError::INVALID_STATE );
171 : 0 : return 0.0;
172 : : }
173 [ + - ]: 130 : patch_set->set_mesh( mesh );
174 : :
175 [ + - ]: 260 : std::vector< PatchSet::PatchHandle > patch_list;
176 [ + - ]: 130 : patch_set->get_patch_handles( patch_list, err );
177 [ + - ][ - + ]: 130 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
178 : :
179 : : // check for inverted elements when using a barrietr target metric
180 : 130 : TQualityMetric* tqm_ptr = NULL;
181 : 130 : TMetricBarrier* tm_ptr = NULL;
182 : 130 : AWMetricBarrier* awm_ptr = NULL;
183 : 130 : ElemSampleQM* sample_qm_ptr = NULL;
184 : 130 : QualityMetric* qm_ptr = NULL;
185 : 130 : TQualityMetric* pmeanp_ptr = NULL;
186 : 130 : ElementMaxQM* elem_max_ptr = NULL;
187 : 130 : ElementAvgQM* elem_avg_ptr = NULL;
188 : 130 : ElementPMeanP* elem_pmeanp_ptr = NULL;
189 : :
190 [ + - ][ + + ]: 130 : ObjectiveFunctionTemplate* of_ptr = dynamic_cast< ObjectiveFunctionTemplate* >( obj_func.get_objective_function() );
191 [ + + ][ + - ]: 130 : if( of_ptr ) qm_ptr = of_ptr->get_quality_metric();
192 [ + + ]: 130 : if( qm_ptr )
193 : : {
194 [ - + ]: 123 : pmeanp_ptr = dynamic_cast< TQualityMetric* >( qm_ptr ); // PMeanP case
195 [ - + ]: 123 : elem_max_ptr = dynamic_cast< ElementMaxQM* >( qm_ptr );
196 [ - + ]: 123 : elem_avg_ptr = dynamic_cast< ElementAvgQM* >( qm_ptr );
197 [ - + ]: 123 : elem_pmeanp_ptr = dynamic_cast< ElementPMeanP* >( qm_ptr );
198 : : }
199 [ + + ]: 130 : if( elem_max_ptr )
200 [ + - ]: 2 : sample_qm_ptr = elem_max_ptr->get_quality_metric();
201 [ + + ]: 128 : else if( elem_pmeanp_ptr )
202 [ + - ]: 34 : sample_qm_ptr = elem_pmeanp_ptr->get_quality_metric();
203 [ - + ]: 94 : else if( elem_avg_ptr )
204 [ # # ]: 0 : sample_qm_ptr = elem_avg_ptr->get_quality_metric();
205 [ + + ]: 94 : else if( pmeanp_ptr )
206 : : {
207 [ + - ][ - + ]: 45 : tm_ptr = dynamic_cast< TMetricBarrier* >( pmeanp_ptr->get_target_metric() );
208 [ + - ][ - + ]: 45 : awm_ptr = dynamic_cast< AWMetricBarrier* >( pmeanp_ptr->get_target_metric() );
209 : : }
210 : :
211 [ + + ][ + + ]: 130 : if( sample_qm_ptr || pmeanp_ptr )
212 : : {
213 [ + + ]: 81 : if( !pmeanp_ptr )
214 : : {
215 [ - + ]: 36 : tqm_ptr = dynamic_cast< TQualityMetric* >( sample_qm_ptr );
216 [ + + ]: 36 : if( tqm_ptr )
217 : : {
218 [ + - ][ - + ]: 35 : tm_ptr = dynamic_cast< TMetricBarrier* >( tqm_ptr->get_target_metric() );
219 [ + - ][ - + ]: 35 : awm_ptr = dynamic_cast< AWMetricBarrier* >( tqm_ptr->get_target_metric() );
220 : : }
221 : : }
222 : : else
223 : : {
224 : 45 : tqm_ptr = dynamic_cast< TQualityMetric* >( pmeanp_ptr );
225 : : }
226 : :
227 [ + + ][ + + ]: 81 : if( tqm_ptr && ( tm_ptr || awm_ptr ) )
[ - + ]
228 : : {
229 : : // check for inverted elements
230 [ + - ]: 31 : this->initialize( patch, err );
231 [ + - ]: 31 : std::vector< size_t > handles;
232 : :
233 : : // set up patch data
234 : 31 : std::vector< PatchSet::PatchHandle >::iterator patch_iter = patch_list.begin();
235 [ + - ][ + + ]: 92 : while( patch_iter != patch_list.end() )
[ + - + ]
236 : : {
237 [ - + ]: 61 : do
238 : : {
239 [ + - ][ + - ]: 61 : patch_set->get_patch( *patch_iter, patch_elements, patch_vertices, err );
240 [ + - ][ - + ]: 61 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
241 [ + - ]: 61 : ++patch_iter;
242 [ - + ][ # # ]: 61 : } while( patch_elements.empty() && patch_iter != patch_list.end() );
[ # # ][ + - ]
[ # # ]
243 : :
244 [ + - ]: 61 : patch.set_mesh_entities( patch_elements, patch_vertices, err );
245 [ + - ][ - + ]: 61 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
246 : :
247 [ + - ]: 61 : qm_ptr->get_evaluations( patch, handles, true, err ); // MSQ_ERRFALSE(err);
248 : :
249 : : // do actual check for inverted elements
250 : 61 : std::vector< size_t >::const_iterator i;
251 : : double tvalue;
252 [ + - ][ + - ]: 33783 : for( i = handles.begin(); i != handles.end(); ++i )
[ + - ][ + + ]
253 : : {
254 [ + - ][ + - ]: 33724 : bool result = tqm_ptr->evaluate( patch, *i, tvalue, err );
255 [ + - ][ + + ]: 33724 : if( MSQ_CHKERR( err ) || !result ) return false; // inverted element detected
[ + - ][ - + ]
[ - + ][ + + ]
256 : : }
257 : 31 : }
258 : : }
259 : : }
260 : :
261 : : // Get termination criteria
262 [ + - ]: 128 : outer_crit = this->get_outer_termination_criterion();
263 [ + - ]: 128 : inner_crit = this->get_inner_termination_criterion();
264 [ - + ]: 128 : if( outer_crit == 0 )
265 : : {
266 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Termination Criterion pointer is Null", MsqError::INVALID_STATE );
267 : 0 : return 0.;
268 : : }
269 [ - + ]: 128 : if( inner_crit == 0 )
270 : : {
271 : : MSQ_SETERR( err )
272 [ # # ][ # # ]: 0 : ( "Termination Criterion pointer for inner loop is Null", MsqError::INVALID_STATE );
273 : 0 : return 0.;
274 : : }
275 : :
276 : : // Set Termination Criterion defaults if no other Criterion is set
277 [ + - ][ + + ]: 128 : if( !outer_crit->criterion_is_set() ) outer_crit->add_iteration_limit( 1 );
[ + - ]
278 [ + - ][ + + ]: 128 : if( !inner_crit->criterion_is_set() ) inner_crit->add_iteration_limit( 10 );
[ + - ]
279 : :
280 : : // If using a local patch, suppress output of inner termination criterion
281 [ + + ]: 128 : if( patch_list.size() > 1 )
282 [ + - ]: 34 : inner_crit->set_debug_output_level( 3 );
283 : : else
284 : 94 : one_patch = true;
285 : :
286 [ + + ]: 128 : if( jacobiOpt )
287 : : {
288 [ + - ]: 1 : coord_tag = get_jacobi_coord_tag( mesh, err );
289 [ + - ][ - + ]: 1 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
290 : 1 : coord_tag_ptr = &coord_tag;
291 : : }
292 : :
293 : : // Initialize outer loop
294 : :
295 [ + - ]: 128 : this->initialize( patch, err );
296 [ + - ][ - + ]: 128 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
297 : :
298 [ + - ]: 128 : valid = obj_func.initialize( mesh_and_domain, settings, patch_set, err );
299 [ + - ][ - + ]: 128 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
300 [ - + ]: 128 : if( !valid )
301 : : {
302 : : MSQ_SETERR( err )
303 : : ( "ObjectiveFunction initialization failed. Mesh "
304 : : "invalid at one or more sample points.",
305 [ # # ][ # # ]: 0 : MsqError::INVALID_MESH );
306 : 0 : goto ERROR;
307 : : }
308 : :
309 [ + - ]: 128 : outer_crit->reset_outer( mesh, domain, obj_func, settings, err );
310 [ + - ][ - + ]: 128 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
311 : :
312 : : // if only one patch, get the patch now
313 [ + + ]: 128 : if( one_patch )
314 : : {
315 [ + - ][ + - ]: 94 : patch_set->get_patch( patch_list[0], patch_elements, patch_vertices, err );
316 [ + - ][ - + ]: 94 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
317 [ + - ]: 94 : patch.set_mesh_entities( patch_elements, patch_vertices, err );
318 [ + - ][ - + ]: 94 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
319 : : }
320 : :
321 : : // Loop until outer termination criterion is met
322 : 128 : inner_crit_terminated = false;
323 [ + - ][ + + ]: 957 : while( !outer_crit->terminate() )
324 : : {
325 [ - + ]: 832 : if( inner_crit_terminated )
326 : : {
327 : : MSQ_SETERR( err )
328 : : ( "Inner termiation criterion satisfied for all patches "
329 : : "without meeting outer termination criterion. This is "
330 : : "an infinite loop. Aborting.",
331 [ # # ][ # # ]: 0 : MsqError::INVALID_STATE );
332 : 3 : break;
333 : : }
334 : 832 : inner_crit_terminated = true;
335 : 832 : all_culled = true;
336 : :
337 : 832 : int num_patches = 0;
338 : : // Loop over each patch
339 : 832 : std::vector< PatchSet::PatchHandle >::iterator p_iter = patch_list.begin();
340 [ + - ][ + + ]: 830780 : while( p_iter != patch_list.end() )
341 : : {
342 [ + + ]: 830141 : if( !one_patch )
343 : : { // if only one patch (global) re-use the previous one
344 : : // loop until we get a non-empty patch. patch will be empty
345 : : // for culled vertices with element-on-vertex patches
346 [ + + ]: 951183 : do
347 : : {
348 [ + - ][ + - ]: 890463 : patch_set->get_patch( *p_iter, patch_elements, patch_vertices, err );
349 [ + - ][ - + ]: 890463 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
350 [ + - ]: 890463 : ++p_iter;
351 [ + + ][ + - ]: 951183 : } while( patch_elements.empty() && p_iter != patch_list.end() );
[ + + ][ + + ]
[ # # ]
352 : :
353 [ + + ]: 829936 : if( patch_elements.empty() )
354 : : { // no more non-culled vertices
355 : : if( 0 ) std::cout << "P[" << get_parallel_rank() << "] tmp srk all vertices culled." << std::endl;
356 : 193 : break;
357 : : }
358 : :
359 : 829743 : all_culled = false;
360 [ + - ]: 829743 : patch.set_mesh_entities( patch_elements, patch_vertices, err );
361 [ + - ][ - + ]: 829743 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
362 : : }
363 : : else
364 : : {
365 [ + - ]: 205 : ++p_iter;
366 : 205 : all_culled = false;
367 : : }
368 : :
369 : 829948 : ++num_patches;
370 : :
371 : : // Initialize for inner iteration
372 : :
373 [ + - ]: 829948 : this->initialize_mesh_iteration( patch, err );
374 [ + - ][ - + ]: 829948 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
375 : :
376 [ + - ]: 829948 : obj_func.reset();
377 : :
378 [ + - ]: 829948 : outer_crit->reset_patch( patch, err );
379 [ + - ][ - + ]: 829948 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
380 : :
381 [ + - ]: 829948 : inner_crit->reset_inner( patch, obj_func, err );
382 [ + - ][ - + ]: 829948 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
383 : :
384 [ + - ]: 829948 : inner_crit->reset_patch( patch, err );
385 [ + - ][ - + ]: 829948 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
386 : :
387 : : // Don't even call optimizer if inner termination
388 : : // criterion has already been met.
389 [ + - ][ + + ]: 829948 : if( !inner_crit->terminate() )
390 : : {
391 : 829946 : inner_crit_terminated = false;
392 : :
393 : : // Call optimizer - should loop on inner_crit->terminate()
394 [ + - ]: 829946 : this->optimize_vertex_positions( patch, err );
395 [ + - ][ - + ]: 829946 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
396 : :
397 : : // Update for changes during inner iteration
398 : : // (during optimizer loop)
399 : :
400 [ + - ]: 829946 : outer_crit->accumulate_patch( patch, err );
401 [ + - ][ - + ]: 829946 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
402 : :
403 [ + - ]: 829946 : inner_crit->cull_vertices( patch, obj_func, err );
404 [ + - ][ - + ]: 829946 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
405 : :
406 : : // FIXME
407 : : if( 0 )
408 : : {
409 : : inner_crit->cull_vertices_global( patch, mesh, domain, settings, obj_func, err );
410 : : if( MSQ_CHKERR( err ) ) goto ERROR;
411 : : }
412 : :
413 [ + - ]: 829946 : patch.update_mesh( err, coord_tag_ptr );
414 [ + - ][ - + ]: 829946 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
415 : : }
416 : : }
417 : :
418 [ + + ][ + - ]: 832 : if( jacobiOpt ) commit_jacobi_coords( coord_tag, mesh, err );
419 : :
420 [ + - ]: 832 : this->terminate_mesh_iteration( patch, err );
421 [ + - ][ - + ]: 832 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
422 : :
423 [ + - ]: 832 : outer_crit->accumulate_outer( mesh, domain, obj_func, settings, err );
424 [ + - ][ - + ]: 832 : if( MSQ_CHKERR( err ) ) goto ERROR;
[ # # ][ # # ]
[ - + ]
425 : :
426 [ + + ]: 832 : if( all_culled ) break;
427 : : }
428 : :
429 : : ERROR:
430 [ + + ][ + - ]: 128 : if( jacobiOpt ) mesh->tag_delete( coord_tag, err );
431 : :
432 : : // call the criteria's cleanup funtions.
433 [ + - ][ + - ]: 128 : if( outer_crit ) outer_crit->cleanup( mesh, domain, err );
434 [ + - ][ + - ]: 128 : if( inner_crit ) inner_crit->cleanup( mesh, domain, err );
435 : : // call the optimization cleanup function.
436 [ + - ]: 128 : this->cleanup();
437 : :
438 : 258 : return 0.;
439 : : }
440 : :
441 : 0 : static void checkpoint_bytes( Mesh* mesh, std::vector< unsigned char >& saved_bytes, MsqError& err )
442 : : {
443 [ # # ]: 0 : std::vector< Mesh::VertexHandle > vertexHandlesArray;
444 [ # # ][ # # ]: 0 : mesh->get_all_vertices( vertexHandlesArray, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
445 [ # # ]: 0 : saved_bytes.resize( vertexHandlesArray.size() );
446 [ # # ][ # # ]: 0 : mesh->vertices_get_byte( &vertexHandlesArray[0], &saved_bytes[0], vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
447 : : }
448 : :
449 : 0 : static void restore_bytes( Mesh* mesh, std::vector< unsigned char >& saved_bytes, MsqError& err )
450 : : {
451 [ # # ]: 0 : std::vector< Mesh::VertexHandle > vertexHandlesArray;
452 [ # # ][ # # ]: 0 : mesh->get_all_vertices( vertexHandlesArray, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
453 [ # # ][ # # ]: 0 : mesh->vertices_set_byte( &vertexHandlesArray[0], &saved_bytes[0], vertexHandlesArray.size(), err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
454 : : }
455 : :
456 : 0 : static void save_or_restore_debug_state( bool save )
457 : : {
458 : : static bool debug[3] = { false, false, false };
459 [ # # ]: 0 : if( save )
460 : : {
461 : 0 : debug[0] = MsqDebug::get( 1 );
462 : 0 : debug[1] = MsqDebug::get( 2 );
463 : 0 : debug[2] = MsqDebug::get( 3 );
464 : : }
465 : : else
466 : : {
467 [ # # ]: 0 : if( debug[0] ) MsqDebug::enable( 1 );
468 [ # # ]: 0 : if( debug[1] ) MsqDebug::enable( 2 );
469 [ # # ]: 0 : if( debug[2] ) MsqDebug::enable( 3 );
470 : : }
471 : 0 : }
472 : :
473 : : /*! \brief Improves the quality of the MeshSet, calling some
474 : : methods specified in a class derived from VertexMover
475 : :
476 : : \param const MeshSet &: this MeshSet is looped over. Only the
477 : : mutable data members are changed (such as currentVertexInd).
478 : : */
479 : 0 : double VertexMover::loop_over_mesh( ParallelMesh* mesh, MeshDomain* domain, const Settings* settings, MsqError& err )
480 : : {
481 [ # # ]: 0 : std::vector< size_t > junk;
482 : : Mesh::VertexHandle vertex_handle;
483 : 0 : TagHandle coord_tag = 0; // store uncommitted coords for jacobi optimization
484 : 0 : TagHandle* coord_tag_ptr = 0;
485 : 0 : int outer_iter = 0;
486 : 0 : int inner_iter = 0;
487 : 0 : bool one_patch = false;
488 : :
489 : : // Clear culling flag, set hard fixed flag, etc on all vertices
490 [ # # ][ # # ]: 0 : MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( (Mesh*)mesh, 0 );
491 [ # # ]: 0 : initialize_vertex_byte( &mesh_and_domain, settings, err );
492 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
493 : :
494 : : // Get the patch data to use for the first iteration
495 [ # # ]: 0 : OFEvaluator& obj_func = get_objective_function_evaluator();
496 : :
497 [ # # ]: 0 : PatchData patch;
498 [ # # ][ # # ]: 0 : patch.set_mesh( (Mesh*)mesh );
499 [ # # ]: 0 : patch.set_domain( domain );
500 [ # # ]: 0 : patch.attach_settings( settings );
501 : :
502 [ # # ]: 0 : ParallelHelper* helper = mesh->get_parallel_helper();
503 [ # # ]: 0 : if( !helper )
504 : : {
505 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No ParallelHelper instance", MsqError::INVALID_STATE );
506 : 0 : return 0;
507 : : }
508 : :
509 [ # # ]: 0 : helper->smoothing_init( err );
510 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
511 : :
512 : : bool inner_crit_terminated, all_culled;
513 [ # # ]: 0 : std::vector< Mesh::VertexHandle > patch_vertices;
514 [ # # ]: 0 : std::vector< Mesh::ElementHandle > patch_elements;
515 [ # # ]: 0 : std::vector< Mesh::VertexHandle > fixed_vertices;
516 [ # # ]: 0 : std::vector< Mesh::VertexHandle > free_vertices;
517 : :
518 : : // Get termination criteria
519 [ # # ]: 0 : TerminationCriterion* outer_crit = this->get_outer_termination_criterion();
520 [ # # ]: 0 : TerminationCriterion* inner_crit = this->get_inner_termination_criterion();
521 [ # # ]: 0 : if( outer_crit == 0 )
522 : : {
523 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Termination Criterion pointer is Null", MsqError::INVALID_STATE );
524 : 0 : return 0.;
525 : : }
526 [ # # ]: 0 : if( inner_crit == 0 )
527 : : {
528 : : MSQ_SETERR( err )
529 [ # # ][ # # ]: 0 : ( "Termination Criterion pointer for inner loop is Null", MsqError::INVALID_STATE );
530 : 0 : return 0.;
531 : : }
532 : :
533 [ # # ]: 0 : PatchSet* patch_set = get_patch_set();
534 [ # # ]: 0 : if( !patch_set )
535 : : {
536 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No PatchSet for QualityImprover!", MsqError::INVALID_STATE );
537 : 0 : return 0.0;
538 : : }
539 [ # # ][ # # ]: 0 : patch_set->set_mesh( (Mesh*)mesh );
540 : :
541 [ # # ]: 0 : std::vector< PatchSet::PatchHandle > patch_list;
542 [ # # ]: 0 : patch_set->get_patch_handles( patch_list, err );
543 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
544 : :
545 [ # # ]: 0 : if( patch_list.size() > 1 )
546 [ # # ]: 0 : inner_crit->set_debug_output_level( 3 );
547 : : else
548 : 0 : one_patch = true;
549 : :
550 [ # # ]: 0 : if( jacobiOpt )
551 : : {
552 [ # # ][ # # ]: 0 : coord_tag = get_jacobi_coord_tag( mesh, err );
553 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
554 : 0 : coord_tag_ptr = &coord_tag;
555 : : }
556 : :
557 : : // parallel error checking
558 [ # # ]: 0 : MsqError perr;
559 : 0 : bool pdone = false;
560 : :
561 : : #define PERROR_COND continue
562 : :
563 : : // Initialize outer loop
564 : :
565 [ # # ]: 0 : this->initialize( patch, err );
566 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) ) { MSQ_SETERR( perr )( "initialize patch", MsqError::INVALID_STATE ); } // goto ERROR;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
567 : :
568 [ # # ][ # # ]: 0 : MeshDomainAssoc mesh_and_domain2 = MeshDomainAssoc( (Mesh*)mesh, domain );
569 [ # # ]: 0 : obj_func.initialize( &mesh_and_domain2, settings, patch_set, err );
570 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) ) { MSQ_SETERR( perr )( "initialize obj_func", MsqError::INVALID_STATE ); } // goto ERROR;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
571 : :
572 [ # # ][ # # ]: 0 : outer_crit->reset_outer( (Mesh*)mesh, domain, obj_func, settings, err );
573 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) ) { MSQ_SETERR( perr )( "reset_outer", MsqError::INVALID_STATE ); } // goto ERROR;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
574 : :
575 : : // Loop until outer termination criterion is met
576 : 0 : inner_crit_terminated = false;
577 : 0 : all_culled = false;
578 : : for( ;; )
579 : : {
580 : : if( 0 )
581 : : std::cout << "P[" << get_parallel_rank() << "] tmp srk inner_iter= " << inner_iter
582 : : << " outer_iter= " << outer_iter << std::endl;
583 : :
584 : : // PERROR:
585 : :
586 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( perr ) )
[ # # ][ # # ]
[ # # ]
587 : : {
588 [ # # ][ # # ]: 0 : std::cout << "P[" << get_parallel_rank() << "] VertexMover::loop_over_mesh found parallel error: " << perr
[ # # ][ # # ]
[ # # ]
589 [ # # ][ # # ]: 0 : << "\n quitting... pdone= " << pdone << std::endl;
[ # # ]
590 : 0 : pdone = true;
591 : : }
592 : :
593 [ # # ]: 0 : helper->communicate_any_true( pdone, err );
594 : :
595 : : if( 0 )
596 : : std::cout << "P[" << get_parallel_rank() << "] tmp srk inner_iter= " << inner_iter
597 : : << " outer_iter= " << outer_iter << " pdone= " << pdone << std::endl;
598 : :
599 [ # # ]: 0 : if( pdone )
600 : : {
601 [ # # ][ # # ]: 0 : std::cout << "P[" << get_parallel_rank()
[ # # ]
602 [ # # ][ # # ]: 0 : << "] VertexMover::loop_over_mesh found parallel error, quitting... pdone= " << pdone
603 [ # # ]: 0 : << std::endl;
604 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "PARALLEL ERROR", MsqError::PARALLEL_ERROR );
605 : 0 : break;
606 : : }
607 : :
608 : 0 : ++outer_iter;
609 : :
610 : : /// [email protected] 1/19/12: the logic here was changed so that all proc's must agree
611 : : /// on the values used for outer and inner termination before the iteration is stopped.
612 : : /// Previously, the ParallelHelper::communicate_all_true method returned true if any
613 : : /// proc sent it a true value, which seems to be a bug, at least in the name of the
614 : : /// method. The method has been changed to return true only if all proc's values are true.
615 : : /// In the previous version, this meant that if any proc hit its inner or outer
616 : : /// termination criterion, the loop was exited, and thus some parts of the mesh
617 : : /// are potentially left unconverged. Also, if the outer criterion was satisfied on
618 : : /// part of the mesh (say a uniform part), the iterations were not executed at all.
619 : : /// Also, changed name of "did_some" to "inner_crit_terminated", and flipped its boolean
620 : : /// value to be consistent with the name - for readability and for correctness since
621 : : /// we want to communicate a true value through the helper.
622 : :
623 [ # # ]: 0 : bool outer_crit_terminated = outer_crit->terminate();
624 : 0 : bool outer_crit_terminated_local = outer_crit_terminated;
625 [ # # ]: 0 : helper->communicate_all_true( outer_crit_terminated, err );
626 : :
627 : 0 : bool inner_crit_terminated_local = inner_crit_terminated;
628 [ # # ]: 0 : helper->communicate_all_true( inner_crit_terminated, err );
629 : :
630 : 0 : bool all_culled_local = all_culled;
631 [ # # ]: 0 : helper->communicate_all_true( all_culled, err );
632 : :
633 [ # # ][ # # ]: 0 : bool done = all_culled || outer_crit_terminated;
634 [ # # ]: 0 : if( inner_crit_terminated )
635 : : {
636 : : MSQ_SETERR( err )
637 : : ( "Inner termination criterion satisfied for all patches "
638 : : "without meeting outer termination criterion. This is "
639 : : "an infinite loop. Aborting.",
640 [ # # ][ # # ]: 0 : MsqError::INVALID_STATE );
641 : 0 : done = true;
642 [ # # ]: 0 : helper->communicate_any_true( done, err );
643 : : }
644 : :
645 : 0 : bool local_done = done;
646 : :
647 [ # # ]: 0 : helper->communicate_all_true( done, err );
648 : :
649 : : if( 0 )
650 : : std::cout << "P[" << get_parallel_rank() << "] tmp srk done= " << done << " local_done= " << local_done
651 : : << " all_culled= " << all_culled << " outer_crit->terminate()= " << outer_crit->terminate()
652 : : << " outer_term= " << outer_crit_terminated
653 : : << " outer_term_local= " << outer_crit_terminated_local
654 : : << " inner_crit_terminated = " << inner_crit_terminated
655 : : << " inner_crit_terminated_local = " << inner_crit_terminated_local
656 : : << " all_culled = " << all_culled << " all_culled_local = " << all_culled_local << std::endl;
657 : :
658 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) ) { MSQ_SETERR( perr )( "loop start", MsqError::INVALID_STATE ); } // goto ERROR;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
659 : :
660 [ # # ]: 0 : if( done ) break;
661 : :
662 : 0 : inner_crit_terminated = true;
663 : 0 : all_culled = true;
664 : :
665 : : ///*** smooth the interior ***////
666 : :
667 : : // get the fixed vertices (i.e. the ones *not* part of the first independent set)
668 [ # # ]: 0 : helper->compute_first_independent_set( fixed_vertices );
669 : :
670 : : // sort the fixed vertices
671 [ # # ]: 0 : std::sort( fixed_vertices.begin(), fixed_vertices.end() );
672 : :
673 : : // Loop over each patch
674 : : if( 0 && MSQ_DBG( 2 ) )
675 : : std::cout << "P[" << get_parallel_rank() << "] tmp srk number of patches = " << patch_list.size()
676 : : << " inner_iter= " << inner_iter << " outer_iter= " << outer_iter
677 : : << " inner.globalInvertedCount = " << inner_crit->globalInvertedCount
678 : : << " outer.globalInvertedCount = " << outer_crit->globalInvertedCount
679 : : << " inner.patchInvertedCount = " << inner_crit->patchInvertedCount
680 : : << " outer.patchInvertedCount = " << outer_crit->patchInvertedCount << std::endl;
681 : :
682 [ # # ]: 0 : save_or_restore_debug_state( true );
683 : : // MsqDebug::disable_all();
684 : :
685 : 0 : int num_patches = 0;
686 : :
687 : 0 : std::vector< PatchSet::PatchHandle >::iterator p_iter = patch_list.begin();
688 [ # # ][ # # ]: 0 : while( p_iter != patch_list.end() )
689 : : {
690 : :
691 : : // loop until we get a non-empty patch. patch will be empty
692 : : // for culled vertices with element-on-vertex patches
693 [ # # ]: 0 : do
694 : : {
695 [ # # ][ # # ]: 0 : patch_set->get_patch( *p_iter, patch_elements, patch_vertices, err );
696 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
697 : : {
698 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "get_patch", MsqError::INVALID_STATE );
699 : 0 : PERROR_COND;
700 : : } // goto ERROR;
701 [ # # ]: 0 : ++p_iter;
702 [ # # ][ # # ]: 0 : } while( patch_elements.empty() && p_iter != patch_list.end() );
[ # # ][ # # ]
[ # # ]
703 : :
704 [ # # ]: 0 : if( patch_elements.empty() )
705 : : { // no more non-culled vertices
706 : : if( 0 ) std::cout << "P[" << get_parallel_rank() << "] tmp srk all vertices culled." << std::endl;
707 : 0 : break;
708 : : }
709 : :
710 [ # # ]: 0 : if( patch_vertices.empty() ) // global patch hack (means all mesh vertices)
711 [ # # ]: 0 : { mesh->get_all_vertices( patch_vertices, err ); }
712 : :
713 : 0 : free_vertices.clear();
714 : :
715 [ # # ]: 0 : for( size_t i = 0; i < patch_vertices.size(); ++i )
716 [ # # ][ # # ]: 0 : if( !std::binary_search( fixed_vertices.begin(), fixed_vertices.end(), patch_vertices[i] ) )
[ # # ]
717 [ # # ][ # # ]: 0 : free_vertices.push_back( patch_vertices[i] );
718 : :
719 [ # # ]: 0 : if( free_vertices.empty() )
720 : : { // all vertices were fixed -> skip patch
721 : 0 : continue;
722 : : }
723 : :
724 : 0 : ++num_patches;
725 : 0 : all_culled = false;
726 [ # # ]: 0 : patch.set_mesh_entities( patch_elements, free_vertices, err );
727 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
728 : : {
729 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "set_mesh_entities", MsqError::INVALID_STATE );
730 : 0 : PERROR_COND;
731 : : } // goto ERROR;
732 : :
733 : : // Initialize for inner iteration
734 : :
735 [ # # ]: 0 : this->initialize_mesh_iteration( patch, err );
736 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
737 : : {
738 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "initialize_mesh_iteration", MsqError::INVALID_STATE );
739 : 0 : PERROR_COND;
740 : : } // goto ERROR;
741 : :
742 [ # # ]: 0 : obj_func.reset();
743 : :
744 [ # # ]: 0 : outer_crit->reset_patch( patch, err );
745 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
746 : : {
747 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "reset_patch outer", MsqError::INVALID_STATE );
748 : 0 : PERROR_COND;
749 : : } // goto ERROR;
750 : :
751 [ # # ]: 0 : inner_crit->reset_inner( patch, obj_func, err );
752 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
753 : : {
754 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "reset_inner", MsqError::INVALID_STATE );
755 : 0 : PERROR_COND;
756 : : } // goto ERROR;
757 : :
758 [ # # ]: 0 : inner_crit->reset_patch( patch, err );
759 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
760 : : {
761 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "inner reset_patch", MsqError::INVALID_STATE );
762 : 0 : PERROR_COND;
763 : : } // goto ERROR;
764 : :
765 : : // Don't even call optimizer if inner termination
766 : : // criterion has already been met.
767 [ # # ][ # # ]: 0 : if( !inner_crit->terminate() )
768 : : {
769 : 0 : inner_crit_terminated = false;
770 [ # # ]: 0 : if( one_patch ) ++inner_iter;
771 : :
772 : : // Call optimizer - should loop on inner_crit->terminate()
773 : : // size_t num_vert=patch.num_free_vertices();
774 : : // std::cout << "P[" << get_parallel_rank() << "] tmp srk VertexMover num_vert= " <<
775 : : // num_vert << std::endl;
776 : :
777 [ # # ]: 0 : this->optimize_vertex_positions( patch, err );
778 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
779 : : {
780 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "optimize_vertex_positions", MsqError::INVALID_STATE );
781 : 0 : PERROR_COND;
782 : : } // goto ERROR;
783 : :
784 : : // Update for changes during inner iteration
785 : : // (during optimizer loop)
786 : :
787 [ # # ]: 0 : outer_crit->accumulate_patch( patch, err );
788 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
789 : : {
790 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "outer accumulate_patch", MsqError::INVALID_STATE );
791 : 0 : PERROR_COND;
792 : : } // goto ERROR;
793 : :
794 [ # # ]: 0 : inner_crit->cull_vertices( patch, obj_func, err );
795 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
796 : : {
797 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "inner cull_vertices", MsqError::INVALID_STATE );
798 : 0 : PERROR_COND;
799 : : } // goto ERROR;
800 : :
801 : : // experimental...
802 : : if( 0 )
803 : : {
804 : : inner_crit->cull_vertices_global( patch, mesh, domain, settings, obj_func, err );
805 : : if( MSQ_CHKERR( err ) )
806 : : {
807 : : MSQ_SETERR( perr )( "cull_vertices_global", MsqError::INVALID_STATE );
808 : : PERROR_COND;
809 : : } // goto ERROR;
810 : : }
811 : :
812 [ # # ]: 0 : patch.update_mesh( err, coord_tag_ptr );
813 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
814 : : {
815 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "update_mesh", MsqError::INVALID_STATE );
816 : 0 : PERROR_COND;
817 : : } // goto ERROR;
818 : : }
819 : : } // while(p_iter....
820 : :
821 [ # # ]: 0 : save_or_restore_debug_state( false );
822 : :
823 : : if( 1 )
824 : : {
825 : 0 : bool pdone_inner = false;
826 : :
827 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( perr ) )
[ # # ][ # # ]
[ # # ]
828 : : {
829 [ # # ][ # # ]: 0 : std::cout << "P[" << get_parallel_rank()
[ # # ]
830 [ # # ][ # # ]: 0 : << "] VertexMover::loop_over_mesh found parallel error: " << perr
831 [ # # ][ # # ]: 0 : << "\n quitting... pdone_inner= " << pdone_inner << std::endl;
[ # # ]
832 : 0 : pdone_inner = true;
833 : : }
834 : :
835 [ # # ]: 0 : helper->communicate_any_true( pdone_inner, err );
836 : :
837 : : if( 0 )
838 : : std::cout << "P[" << get_parallel_rank() << "] tmp srk inner_iter= " << inner_iter
839 : : << " outer_iter= " << outer_iter << " pdone_inner= " << pdone_inner << std::endl;
840 : :
841 [ # # ]: 0 : if( pdone_inner )
842 : : {
843 : : if( 0 )
844 : : std::cout << "P[" << get_parallel_rank()
845 : : << "] tmp srk found parallel error, quitting... pdone_inner= " << pdone_inner
846 : : << std::endl;
847 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "PARALLEL ERROR", MsqError::PARALLEL_ERROR );
848 : 0 : break;
849 : : }
850 : : }
851 : :
852 : : /// [email protected] save vertex bytes since boundary smoothing changes them
853 [ # # ]: 0 : std::vector< unsigned char > saved_bytes;
854 [ # # ][ # # ]: 0 : checkpoint_bytes( mesh, saved_bytes, err );
855 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
856 : : {
857 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "checkpoint_bytes ", MsqError::INVALID_STATE );
858 : 0 : PERROR_COND;
859 : : } // goto ERROR;
860 : :
861 [ # # ]: 0 : helper->communicate_first_independent_set( err );
862 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
863 : : {
864 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "communicate_first_independent_set ", MsqError::INVALID_STATE );
865 : 0 : PERROR_COND;
866 : : } // goto ERROR;
867 : :
868 : : ///*** smooth the boundary ***////
869 [ # # ]: 0 : save_or_restore_debug_state( true );
870 [ # # ]: 0 : MsqDebug::disable_all();
871 : :
872 [ # # ][ # # ]: 0 : while( helper->compute_next_independent_set() )
873 : : {
874 : : // Loop over all boundary elements
875 [ # # ][ # # ]: 0 : while( helper->get_next_partition_boundary_vertex( vertex_handle ) )
876 : : {
877 : :
878 : 0 : patch_vertices.clear();
879 [ # # ]: 0 : patch_vertices.push_back( vertex_handle );
880 : 0 : patch_elements.clear();
881 [ # # ]: 0 : mesh->vertices_get_attached_elements( &vertex_handle, 1, patch_elements, junk, err );
882 : :
883 : 0 : all_culled = false;
884 [ # # ]: 0 : patch.set_mesh_entities( patch_elements, patch_vertices, err );
885 : :
886 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
887 : : {
888 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "set_mesh_entities 2 ", MsqError::INVALID_STATE );
889 : 0 : PERROR_COND;
890 : : } // goto ERROR;
891 : :
892 : : // Initialize for inner iteration
893 : :
894 [ # # ]: 0 : this->initialize_mesh_iteration( patch, err );
895 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
896 : : {
897 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " initialize_mesh_iteration 2", MsqError::INVALID_STATE );
898 : 0 : PERROR_COND;
899 : : } // goto ERROR;
900 : :
901 [ # # ]: 0 : obj_func.reset();
902 : :
903 [ # # ]: 0 : outer_crit->reset_patch( patch, err );
904 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
905 : : {
906 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( "outer reset_patch 2 ", MsqError::INVALID_STATE );
907 : 0 : PERROR_COND;
908 : : } // goto ERROR;
909 : :
910 [ # # ]: 0 : inner_crit->reset_inner( patch, obj_func, err );
911 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
912 : : {
913 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " inner reset_inner 2", MsqError::INVALID_STATE );
914 : 0 : PERROR_COND;
915 : : } // goto ERROR;
916 : :
917 [ # # ]: 0 : inner_crit->reset_patch( patch, err );
918 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
919 : : {
920 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " inner_crit reset_patch 2", MsqError::INVALID_STATE );
921 : 0 : PERROR_COND;
922 : : } // goto ERROR;
923 : :
924 : : // Don't even call optimizer if inner termination
925 : : // criterion has already been met.
926 [ # # ][ # # ]: 0 : if( !inner_crit->terminate() )
927 : : {
928 : 0 : inner_crit_terminated = false;
929 : :
930 : : // Call optimizer - should loop on inner_crit->terminate()
931 [ # # ]: 0 : this->optimize_vertex_positions( patch, err );
932 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
933 : : {
934 : : MSQ_SETERR( perr )
935 [ # # ][ # # ]: 0 : ( " optimize_vertex_positions 2", MsqError::INVALID_STATE );
936 : 0 : PERROR_COND;
937 : : } // goto ERROR;
938 : :
939 : : // Update for changes during inner iteration
940 : : // (during optimizer loop)
941 : :
942 [ # # ]: 0 : outer_crit->accumulate_patch( patch, err );
943 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
944 : : {
945 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " outer accumulate_patch 2", MsqError::INVALID_STATE );
946 : 0 : PERROR_COND;
947 : : } // goto ERROR;
948 : :
949 [ # # ]: 0 : inner_crit->cull_vertices( patch, obj_func, err );
950 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
951 : : {
952 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " inner cull_vertices", MsqError::INVALID_STATE );
953 : 0 : PERROR_COND;
954 : : } // goto ERROR;
955 : :
956 : : // FIXME
957 : : if( 0 )
958 : : {
959 : : inner_crit->cull_vertices_global( patch, mesh, domain, settings, obj_func, err );
960 : : if( MSQ_CHKERR( err ) )
961 : : {
962 : : MSQ_SETERR( perr )
963 : : ( " cull_vertices_global 2 ", MsqError::INVALID_STATE );
964 : : PERROR_COND;
965 : : } // goto ERROR;
966 : : }
967 : :
968 [ # # ]: 0 : patch.update_mesh( err, coord_tag_ptr );
969 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
970 : : {
971 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " update_mesh 2", MsqError::INVALID_STATE );
972 : 0 : PERROR_COND;
973 : : } // goto ERROR;
974 : : }
975 : : }
976 [ # # ]: 0 : helper->communicate_next_independent_set( err );
977 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
978 : : {
979 : : MSQ_SETERR( perr )
980 [ # # ][ # # ]: 0 : ( " communicate_next_independent_set 2", MsqError::INVALID_STATE );
981 : 0 : PERROR_COND;
982 : : } // goto ERROR;
983 : : } // while(helper->compute_next_independent_set())
984 : :
985 [ # # ]: 0 : save_or_restore_debug_state( false );
986 : :
987 : : // if (!get_parallel_rank())
988 : : // std::cout << "P[" << get_parallel_rank() << "] tmp srk num_patches= " << num_patches <<
989 : : // std::endl;
990 : :
991 : : /// [email protected] restore vertex bytes since boundary smoothing changes them
992 [ # # ][ # # ]: 0 : restore_bytes( mesh, saved_bytes, err );
993 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
994 : : {
995 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " restore_bytes", MsqError::INVALID_STATE );
996 : 0 : PERROR_COND;
997 : : } // goto ERROR;
998 : :
999 [ # # ][ # # ]: 0 : if( jacobiOpt ) commit_jacobi_coords( coord_tag, mesh, err );
[ # # ]
1000 : :
1001 [ # # ]: 0 : this->terminate_mesh_iteration( patch, err );
1002 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
1003 : : {
1004 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " terminate_mesh_iteration", MsqError::INVALID_STATE );
1005 : 0 : PERROR_COND;
1006 : : } // goto ERROR;
1007 : :
1008 [ # # ][ # # ]: 0 : outer_crit->accumulate_outer( mesh, domain, obj_func, settings, err );
1009 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
1010 : : {
1011 [ # # ][ # # ]: 0 : MSQ_SETERR( perr )( " outer_crit accumulate_outer", MsqError::INVALID_STATE );
1012 [ # # ]: 0 : PERROR_COND;
1013 : : } // goto ERROR;
1014 : 0 : }
1015 : :
1016 : : // ERROR:
1017 : :
1018 [ # # ][ # # ]: 0 : if( MSQ_CHKERR( err ) )
[ # # ][ # # ]
[ # # ]
1019 : : {
1020 [ # # ][ # # ]: 0 : std::cout << "P[" << get_parallel_rank() << "] VertexMover::loop_over_mesh error = " << err.error_message()
[ # # ][ # # ]
[ # # ][ # # ]
1021 [ # # ]: 0 : << std::endl;
1022 : : }
1023 : :
1024 [ # # ][ # # ]: 0 : if( jacobiOpt ) mesh->tag_delete( coord_tag, err );
1025 : :
1026 : : // call the criteria's cleanup funtions.
1027 [ # # ][ # # ]: 0 : outer_crit->cleanup( mesh, domain, err );MSQ_CHKERR( err );
[ # # ][ # # ]
[ # # ][ # # ]
1028 [ # # ][ # # ]: 0 : inner_crit->cleanup( mesh, domain, err );MSQ_CHKERR( err );
[ # # ][ # # ]
[ # # ][ # # ]
1029 : : // call the optimization cleanup function.
1030 [ # # ]: 0 : this->cleanup();
1031 : : // close the helper
1032 [ # # ][ # # ]: 0 : helper->smoothing_close( err );MSQ_CHKERR( err );
[ # # ][ # # ]
[ # # ]
1033 : :
1034 : 0 : return 0.;
1035 : : }
1036 : :
1037 : 130 : void VertexMover::initialize_queue( MeshDomainAssoc* mesh_and_domain, const Settings* settings, MsqError& err )
1038 : : {
1039 [ - + ][ # # ]: 130 : QualityImprover::initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
[ - + ]
1040 [ - + ][ # # ]: 130 : objFuncEval.initialize_queue( mesh_and_domain, settings, err );MSQ_ERRRTN( err );
[ - + ]
1041 : : }
1042 : :
1043 : 1 : TagHandle VertexMover::get_jacobi_coord_tag( Mesh* mesh, MsqError& err )
1044 : : {
1045 : : // Get tag handle
1046 : 1 : const char tagname[] = "msq_jacobi_temp_coords";
1047 [ + - ][ + - ]: 1 : TagHandle tag = mesh->tag_create( tagname, Mesh::DOUBLE, 3, 0, err );
1048 [ + - ][ - + ]: 1 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
1049 : : /* VertexMover will always delete the tag when it is done, so it is probably
1050 : : best not to accept an existing tag
1051 : : if (err.error_code() == TAG_ALREADY_EXISTS) {
1052 : : err.clear();
1053 : : tag = tag_get( tagname, err ); MSQ_ERRZERO(err);
1054 : : std::string name;
1055 : : Mesh::TagType type;
1056 : : unsigned length;
1057 : : mesh->tag_properties( tag, name, type, length, err ); MSQ_ERRZERO(err);
1058 : : if (type != Mesh::DOUBLE || length != 3) {
1059 : : MSQ_SETERR(err)(TAG_ALREADY_EXISTS,
1060 : : "Tag \"%s\" already exists with invalid type",
1061 : : tagname);
1062 : : return 0;
1063 : : }
1064 : : }
1065 : : */
1066 : :
1067 : : // Initialize tag value with initial vertex coordinates so that
1068 : : // vertices that never get moved end up with the correct coords.
1069 [ + - ]: 1 : std::vector< Mesh::VertexHandle > vertices;
1070 [ + - ]: 1 : mesh->get_all_vertices( vertices, err );
1071 [ + - ][ - + ]: 1 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
1072 : : // remove fixed vertices
1073 : : // to avoid really huge arrays (especially since we need to copy
1074 : : // coords out of annoying MsqVertex class to double array), work
1075 : : // in blocks
1076 : 1 : const size_t blocksize = 4096;
1077 [ + - ]: 2 : std::vector< bool > fixed( blocksize );
1078 [ + - ][ + + ]: 4097 : MsqVertex coords[blocksize];
1079 : : double tag_data[3 * blocksize];
1080 [ + + ]: 2 : for( size_t i = 0; i < vertices.size(); i += blocksize )
1081 : : {
1082 [ + - ]: 1 : size_t count = std::min( blocksize, vertices.size() - i );
1083 : : // remove fixed vertices
1084 [ + - ][ + - ]: 1 : mesh->vertices_get_fixed_flag( &vertices[i], fixed, count, err );
1085 [ + - ][ - + ]: 1 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
1086 : 1 : size_t w = 0;
1087 [ + + ]: 17 : for( size_t j = 0; j < count; ++j )
1088 [ + - ][ + + ]: 16 : if( !fixed[j] ) vertices[i + w++] = vertices[i + j];
[ + - ][ + - ]
1089 : 1 : count = w;
1090 : : // store tag data for free vertices
1091 [ + - ][ + - ]: 1 : mesh->vertices_get_coordinates( &vertices[i], coords, count, err );
1092 [ + - ][ - + ]: 1 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
1093 [ + + ]: 5 : for( size_t j = 0; j < count; ++j )
1094 : : {
1095 [ + - ]: 4 : tag_data[3 * j] = coords[j][0];
1096 [ + - ]: 4 : tag_data[3 * j + 1] = coords[j][1];
1097 [ + - ]: 4 : tag_data[3 * j + 2] = coords[j][2];
1098 : : }
1099 [ + - ][ + - ]: 1 : mesh->tag_set_vertex_data( tag, count, &vertices[i], tag_data, err );
1100 [ + - ][ - + ]: 1 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
1101 : : }
1102 : :
1103 : 2 : return tag;
1104 : : }
1105 : :
1106 : 1 : void VertexMover::commit_jacobi_coords( TagHandle tag, Mesh* mesh, MsqError& err )
1107 : : {
1108 [ + - ]: 1 : Vector3D coords;
1109 [ + - ]: 1 : std::vector< Mesh::VertexHandle > vertices;
1110 [ + - ][ + - ]: 1 : mesh->get_all_vertices( vertices, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
1111 [ + - ][ + - ]: 2 : std::vector< bool > fixed( vertices.size() );
1112 [ + - ][ + - ]: 1 : mesh->vertices_get_fixed_flag( &vertices[0], fixed, vertices.size(), err );
1113 [ + + ][ + - ]: 17 : for( size_t i = 0; i < vertices.size(); ++i )
1114 : : {
1115 [ + - ][ + + ]: 16 : if( !fixed[i] )
1116 : : {
1117 [ + - ][ + - ]: 4 : mesh->tag_get_vertex_data( tag, 1, &vertices[i], coords.to_array(), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
1118 [ + - ][ + - ]: 4 : mesh->vertex_set_coordinates( vertices[i], coords, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
1119 : : }
1120 : 1 : }
1121 : : }
1122 : :
1123 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
