 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /** @example LaplacianSmoother.cpp \n
2 * \brief Perform Laplacian relaxation on a mesh and its dual \n
3 * <b>To run</b>: mpiexec -np <np> LaplacianSmoother [filename]\n
4 *
5 * Briefly, Laplacian relaxation is a technique to smooth out a mesh. The centroid of each cell is
6 * computed from its vertex positions, then vertices are placed at the average of their connected
7 * cells' centroids.
8 *
9 * In the parallel algorithm, an extra ghost layer of cells is exchanged. This allows us to compute
10 * the centroids for boundary cells on each processor where they appear; this eliminates the need
11 * for one round of data exchange (for those centroids) between processors. New vertex positions
12 * must be sent from owning processors to processors sharing those vertices. Convergence is
13 * measured as the maximum distance moved by any vertex.
14 *
15 * In this implementation, a fixed number of iterations is performed. The final mesh is output to
16 * 'laplacianfinal.h5m' in the current directory (H5M format must be used since the file is written
17 * in parallel).
18 *
19 * Usage: mpiexec -n 2 valgrind ./LaplacianSmoother -f input/surfrandomtris-64part.h5m -r 2 -p 2 -n
20 * 25
21 */
22
23 #include <iostream>
24 #include <sstream>
25 #include "moab/Core.hpp"
26 #ifdef MOAB_HAVE_MPI
27 #include "moab/ParallelComm.hpp"
28 #include "MBParallelConventions.h"
29 #endif
30 #include "moab/Skinner.hpp"
31 #include "moab/CN.hpp"
32 #include "moab/ProgOptions.hpp"
33 #include "moab/CartVect.hpp"
34 #include "moab/NestedRefine.hpp"
35 #include "moab/VerdictWrapper.hpp"
36 #include "matrix.h"
37
38 using namespace moab;
39 using namespace std;
40
41 #define WRITE_DEBUG_FILES
42
43 #ifdef MESH_DIR
44 string test_file_name = string( MESH_DIR ) + string( "/surfrandomtris-4part.h5m" );
45 #else
46 string test_file_name = string( "input/surfrandomtris-4part.h5m" );
47 #endif
48
49 #define RC MB_CHK_ERR( rval )
50 #define dbgprint( MSG ) \
51 do \
52 { \
53 if( !global_rank ) std::cerr << ( MSG ) << std::endl; \
54 } while( false )
55
56 ErrorCode perform_laplacian_smoothing( Core* mb,
57 Range& cells,
58 Range& verts,
59 int dim,
60 Tag fixed,
61 bool use_hc = false,
62 bool use_acc = false,
63 int acc_method = 1,
64 int num_its = 10,
65 double rel_eps = 1e-5,
66 double alpha = 0.0,
67 double beta = 0.5,
68 int report_its = 1 );
69
70 ErrorCode hcFilter( Core* mb,
71 moab::ParallelComm* pcomm,
72 moab::Range& verts,
73 int dim,
74 Tag fixed,
75 std::vector< double >& verts_o,
76 std::vector< double >& verts_n,
77 double alpha,
78 double beta );
79
80 ErrorCode laplacianFilter( Core* mb,
81 moab::ParallelComm* pcomm,
82 moab::Range& verts,
83 int dim,
84 Tag fixed,
85 std::vector< double >& verts_o,
86 std::vector< double >& verts_n,
87 bool use_updated = true );
88
89 int main( int argc, char** argv )
90 {
91 int num_its = 10;
92 int num_ref = 0;
93 int num_dim = 2;
94 int report_its = 1;
95 int num_degree = 2;
96 bool use_hc = false;
97 bool use_acc = false;
98 int acc_method = 1;
99 double alpha = 0.5, beta = 0.0;
100 double rel_eps = 1e-5;
101 const int nghostrings = 1;
102 ProgOptions opts;
103 ErrorCode rval;
104 std::stringstream sstr;
105 int global_rank;
106
107 MPI_Init( &argc, &argv );
108 MPI_Comm_rank( MPI_COMM_WORLD, &global_rank );
109
110 // Decipher program options from user
111 opts.addOpt< int >( std::string( "niter,n" ),
112 std::string( "Number of Laplacian smoothing iterations (default=10)" ), &num_its );
113 opts.addOpt< double >( std::string( "eps,e" ),
114 std::string( "Tolerance for the Laplacian smoothing error (default=1e-5)" ), &rel_eps );
115 opts.addOpt< double >( std::string( "alpha" ),
116 std::string( "Tolerance for the Laplacian smoothing error (default=0.0)" ), &alpha );
117 opts.addOpt< double >( std::string( "beta" ),
118 std::string( "Tolerance for the Laplacian smoothing error (default=0.5)" ), &beta );
119 opts.addOpt< int >( std::string( "dim,d" ), std::string( "Topological dimension of the mesh (default=2)" ),
120 &num_dim );
121 opts.addOpt< std::string >( std::string( "file,f" ),
122 std::string( "Input mesh file to smoothen (default=surfrandomtris-4part.h5m)" ),
123 &test_file_name );
124 opts.addOpt< int >(
125 std::string( "nrefine,r" ),
126 std::string( "Number of uniform refinements to perform and apply smoothing cycles (default=1)" ), &num_ref );
127 opts.addOpt< int >( std::string( "ndegree,p" ), std::string( "Degree of uniform refinement (default=2)" ),
128 &num_degree );
129 opts.addOpt< void >( std::string( "humphrey,c" ),
130 std::string( "Use Humphrey’s Classes algorithm to reduce shrinkage of "
131 "Laplacian smoother (default=false)" ),
132 &use_hc );
133 opts.addOpt< void >( std::string( "aitken,d" ),
134 std::string( "Use Aitken \\delta^2 acceleration to improve convergence of "
135 "Lapalace smoothing algorithm (default=false)" ),
136 &use_acc );
137 opts.addOpt< int >( std::string( "acc,a" ),
138 std::string( "Type of vector Aitken process to use for acceleration (default=1)" ),
139 &acc_method );
140
141 opts.parseCommandLine( argc, argv );
142
143 // get MOAB and ParallelComm instances
144 Core* mb = new Core;
145 if( NULL == mb ) return 1;
146 int global_size = 1;
147
148 // get the ParallelComm instance
149 ParallelComm* pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
150 global_size = pcomm->size();
151
152 string roptions, woptions;
153 if( global_size > 1 )
154 { // if reading in parallel, need to tell it how
155 sstr.str( "" );
156 sstr << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
157 "PARALLEL_GHOSTS="
158 << num_dim << ".0." << nghostrings << ";DEBUG_IO=0;DEBUG_PIO=0";
159 roptions = sstr.str();
160 woptions = "PARALLEL=WRITE_PART";
161 }
162
163 // read the file
164 moab::EntityHandle fileset, currset;
165 rval = mb->create_meshset( MESHSET_SET, fileset );
166 RC;
167 currset = fileset;
168 rval = mb->load_file( test_file_name.c_str(), &fileset, roptions.c_str() );
169 RC;
170
171 std::vector< EntityHandle > hsets( num_ref + 1, fileset );
172 if( num_ref )
173 {
174 // Perform uniform refinement of the smoothed mesh
175 #ifdef MOAB_HAVE_MPI
176 NestedRefine* uref = new NestedRefine( mb, pcomm, currset );
177 #else
178 NestedRefine* uref = new NestedRefine( mb, 0, currset );
179 #endif
180
181 std::vector< int > num_degrees( num_ref, num_degree );
182 rval = uref->generate_mesh_hierarchy( num_ref, &num_degrees[0], hsets );
183 RC;
184
185 // Now exchange 1 layer of ghost elements, using vertices as bridge
186 // (we could have done this as part of reading process, using the PARALLEL_GHOSTS read
187 // option)
188 rval = uref->exchange_ghosts( hsets, nghostrings );
189 RC;
190
191 delete uref;
192 }
193
194 for( int iref = 0; iref <= num_ref; ++iref )
195 {
196
197 // specify which set we are currently working on
198 currset = hsets[iref];
199
200 // make tag to specify fixed vertices, since it's input to the algorithm; use a default
201 // value of non-fixed so we only need to set the fixed tag for skin vertices
202 Tag fixed;
203 int def_val = 0;
204 rval = mb->tag_get_handle( "fixed", 1, MB_TYPE_INTEGER, fixed, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );
205 RC;
206
207 // get all vertices and cells
208 Range verts, cells, skin_verts;
209 rval = mb->get_entities_by_type( currset, MBVERTEX, verts );
210 RC;
211 rval = mb->get_entities_by_dimension( currset, num_dim, cells );
212 RC;
213 dbgprint( "Found " << verts.size() << " vertices and " << cells.size() << " elements" );
214
215 // get the skin vertices of those cells and mark them as fixed; we don't want to fix the
216 // vertices on a part boundary, but since we exchanged a layer of ghost cells, those
217 // vertices aren't on the skin locally ok to mark non-owned skin vertices too, I won't move
218 // those anyway use MOAB's skinner class to find the skin
219 Skinner skinner( mb );
220 rval = skinner.find_skin( currset, cells, true, skin_verts );
221 RC; // 'true' param indicates we want vertices back, not cells
222
223 std::vector< int > fix_tag( skin_verts.size(), 1 ); // initialized to 1 to indicate fixed
224 rval = mb->tag_set_data( fixed, skin_verts, &fix_tag[0] );
225 RC;
226 // exchange tags on owned verts for fixed points
227 if( global_size > 1 )
228 {
229 #ifdef MOAB_HAVE_MPI
230 rval = pcomm->exchange_tags( fixed, skin_verts );
231 RC;
232 #endif
233 }
234
235 // now perform the Laplacian relaxation
236 rval = perform_laplacian_smoothing( mb, cells, verts, num_dim, fixed, use_hc, use_acc, acc_method, num_its,
237 rel_eps, alpha, beta, report_its );
238 RC;
239
240 // output file, using parallel write
241 sstr.str( "" );
242 sstr << "LaplacianSmoother_" << iref << ".h5m";
243 rval = mb->write_file( sstr.str().c_str(), "H5M", woptions.c_str(), &currset, 1 );
244 RC;
245
246 // delete fixed tag, since we created it here
247 rval = mb->tag_delete( fixed );
248 RC;
249 }
250
251 delete pcomm;
252 // delete MOAB instance
253 delete mb;
254
255 MPI_Finalize();
256 return 0;
257 }
258
259 ErrorCode perform_laplacian_smoothing( Core* mb,
260 Range& cells,
261 Range& verts,
262 int dim,
263 Tag fixed,
264 bool use_hc,
265 bool use_acc,
266 int acc_method,
267 int num_its,
268 double rel_eps,
269 double alpha,
270 double beta,
271 int report_its )
272 {
273 ErrorCode rval;
274 int global_rank = 0, global_size = 1;
275 int nacc = 2; /* nacc_method: 1 = Method 2 from [1], 2 = Method 3 from [1] */
276 std::vector< double > verts_acc1, verts_acc2, verts_acc3;
277 double rat_theta = rel_eps, rat_alpha = rel_eps, rat_alphaprev = rel_eps;
278 #ifdef MOAB_HAVE_MPI
279 const char* woptions = "PARALLEL=WRITE_PART";
280 #else
281 const char* woptions = "";
282 #endif
283 std::vector< int > fix_tag( verts.size() );
284
285 rval = mb->tag_get_data( fixed, verts, &fix_tag[0] );
286 RC;
287
288 #ifdef MOAB_HAVE_MPI
289 ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
290 global_rank = pcomm->rank();
291 global_size = pcomm->size();
292 #endif
293
294 dbgprint( "-- Starting smoothing cycle --" );
295 // perform Laplacian relaxation:
296 // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
297
298 // get all verts coords into tag; don't need to worry about filtering out fixed verts,
299 // we'll just be setting to their fixed coords
300 std::vector< double > verts_o, verts_n;
301 verts_o.resize( 3 * verts.size(), 0.0 );
302 verts_n.resize( 3 * verts.size(), 0.0 );
303 // std::vector<const void*> vdata(1);
304 // vdata[0] = &verts_n[0];
305 // const int vdatasize = verts_n.size();
306 void* vdata = &verts_n[0];
307 rval = mb->get_coords( verts, &verts_o[0] );
308 RC;
309 const int nbytes = sizeof( double ) * verts_o.size();
310
311 Tag errt, vpost;
312 double def_val[3] = { 0.0, 0.0, 0.0 };
313 rval = mb->tag_get_handle( "error", 1, MB_TYPE_DOUBLE, errt, MB_TAG_CREAT | MB_TAG_DENSE, def_val );
314 RC;
315 rval = mb->tag_get_handle( "vpos", 3, MB_TYPE_DOUBLE, vpost, MB_TAG_CREAT | MB_TAG_DENSE, def_val );
316 RC;
317 // rval = mb->tag_set_by_ptr(vpost, verts, vdata); RC;
318
319 if( use_acc )
320 {
321 verts_acc1.resize( verts_o.size(), 0.0 );
322 verts_acc2.resize( verts_o.size(), 0.0 );
323 verts_acc3.resize( verts_o.size(), 0.0 );
324 memcpy( &verts_acc1[0], &verts_o[0], nbytes );
325 memcpy( &verts_acc2[0], &verts_o[0], nbytes );
326 memcpy( &verts_acc3[0], &verts_o[0], nbytes );
327 }
328
329 // Filter verts down to owned ones and get fixed tag for them
330 Range owned_verts, shared_owned_verts;
331 if( global_size > 1 )
332 {
333 #ifdef MOAB_HAVE_MPI
334 rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );MB_CHK_ERR( rval );
335 #endif
336 }
337 else
338 owned_verts = verts;
339
340 #ifdef MOAB_HAVE_MPI
341 // Get shared owned verts, for exchanging tags
342 rval = pcomm->get_shared_entities( -1, shared_owned_verts, 0, false, true );MB_CHK_ERR( rval );
343 // Workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging
344 // tags for all shared entities
345 if( shared_owned_verts.empty() ) shared_owned_verts.insert( *verts.begin() );
346 #endif
347
348 #ifdef WRITE_DEBUG_FILES
349 {
350 // output file, using parallel write
351 std::stringstream sstr;
352 sstr << "LaplacianSmootherIterate_0.h5m";
353 rval = mb->write_file( sstr.str().c_str(), NULL, woptions );
354 RC;
355 }
356 #endif
357
358 bool check_metrics = true;
359 VerdictWrapper vw( mb );
360 vw.set_size( 1. ); // for relative size measures; maybe it should be user input
361
362 double mxdelta = 0.0, global_max = 0.0;
363 // 2. for num_its iterations:
364 for( int nit = 0; nit < num_its; nit++ )
365 {
366
367 mxdelta = 0.0;
368
369 // 2a. Apply Laplacian Smoothing Filter to Mesh
370 if( use_hc )
371 {
372 rval = hcFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n, alpha, beta );
373 RC;
374 }
375 else
376 {
377 rval = laplacianFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n );
378 RC;
379 }
380
381 if( check_metrics )
382 {
383 bool smooth_success = true;
384 std::vector< double > coords;
385 double jac = 0.0;
386 int num_nodes = 0;
387 for( unsigned ic = 0; ic < cells.size(); ++ic )
388 {
389 EntityHandle ecell = cells[ic];
390 EntityType etype = mb->type_from_handle( ecell );
391 Range everts;
392 rval = mb->get_connectivity( &ecell, 1, everts, num_nodes );
393 RC;
394 coords.resize( num_nodes * 3, 0.0 );
395 for( int iv = 0; iv < num_nodes; ++iv )
396 {
397 const int offset = mb->id_from_handle( everts[iv] ) * 3;
398 for( int ii = 0; ii < 3; ++ii )
399 coords[iv * 3 + ii] = verts_n[offset + ii];
400 }
401 rval = vw.quality_measure( ecell, MB_SHAPE, jac, num_nodes, etype, &coords[0] );
402 RC;
403 if( jac < 1e-8 )
404 {
405 dbgprint( "Inverted element " << ic << " with jacobian = " << jac << "." );
406 smooth_success = false;
407 break;
408 }
409 }
410
411 if( !smooth_success )
412 {
413 verts_o.swap( verts_n );
414 dbgprint( "Found element inversions during smoothing. Stopping iterations." );
415 break;
416 }
417 }
418
419 if( use_acc )
420 {
421
422 // if (acc_method < 5 || nit <= 3) {
423 // memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
424 // memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
425 // memcpy( &verts_acc3[0], &verts_n[0], nbytes );
426 // }
427
428 rat_alphaprev = rat_alpha;
429 for( unsigned i = 0; i < verts_n.size(); ++i )
430 {
431 rat_alpha =
432 std::max( rat_alpha,
433 std::abs( ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc2[i] - verts_acc1[i] ) ) /
434 ( ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] ) ) );
435 }
436 rat_theta = std::abs( rat_alpha / rat_alphaprev - 1.0 );
437
438 if( nit > 3 && ( nit % nacc ) && rat_theta < 1.0 )
439 {
440
441 if( acc_method == 1 )
442 { /* Method 2 from ACCELERATION OF VECTOR SEQUENCES:
443 http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
444 double vnorm = 0.0, den, acc_alpha = 0.0, acc_gamma = 0.0;
445 for( unsigned i = 0; i < verts_n.size(); ++i )
446 {
447 den = ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
448 vnorm += den * den;
449 acc_alpha += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - verts_acc2[i] );
450 acc_gamma += ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] );
451 }
452 for( unsigned i = 0; i < verts_n.size(); ++i )
453 {
454 verts_n[i] = verts_acc2[i] + ( acc_gamma * ( verts_acc3[i] - verts_acc2[i] ) -
455 acc_alpha * ( verts_acc2[i] - verts_acc1[i] ) ) /
456 vnorm;
457 }
458 }
459 else if( acc_method == 2 )
460 { /* Method 3 from ACCELERATION OF VECTOR SEQUENCES:
461 http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
462 double vnorm = 0.0, num = 0.0, den = 0.0, mu = 0.0;
463 for( unsigned i = 0; i < verts_n.size(); ++i )
464 {
465 num +=
466 ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
467 den = ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
468 vnorm += den * den;
469 }
470 mu = num / vnorm;
471 for( unsigned i = 0; i < verts_n.size(); ++i )
472 {
473 verts_n[i] = verts_acc3[i] + mu * ( verts_acc2[i] - verts_acc3[i] );
474 }
475 }
476 else if( acc_method == 3 )
477 { /* Method 5 from ACCELERATION OF VECTOR SEQUENCES:
478 http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
479 double num = 0.0, den = 0.0, mu = 0.0;
480 for( unsigned i = 0; i < verts_n.size(); ++i )
481 {
482 num += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc2[i] - verts_acc1[i] );
483 den +=
484 ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
485 }
486 mu = num / den;
487 for( unsigned i = 0; i < verts_n.size(); ++i )
488 {
489 verts_n[i] = verts_acc3[i] - mu * ( verts_acc3[i] - verts_acc2[i] );
490 }
491 }
492 else if( acc_method == 4 )
493 { /* Method 8 from ACCELERATION OF VECTOR SEQUENCES:
494 http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
495 double num = 0.0, den = 0.0, lambda = 0.0;
496 for( unsigned i = 0; i < verts_n.size(); ++i )
497 {
498 num += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - verts_acc2[i] );
499 den += ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] );
500 }
501 lambda = std::sqrt( num / den );
502 for( unsigned i = 0; i < verts_n.size(); ++i )
503 {
504 verts_n[i] = verts_acc3[i] - lambda / ( lambda - 1.0 ) * ( verts_acc3[i] - verts_acc2[i] );
505 }
506 }
507 else if( acc_method == 5 )
508 { /* Minimum polynomial extrapolation of vector sequenes :
509 https://en.wikipedia.org/wiki/Minimum_polynomial_extrapolation */
510 /* Pseudo-code
511 U=x(:,2:end-1)-x(:,1:end-2);
512 c=-pinv(U)*(x(:,end)-x(:,end-1));
513 c(end+1,1)=1;
514 s=(x(:,2:end)*c)/sum(c);
515 */
516 Matrix U( verts_n.size(), 2 );
517 Vector res( verts_n.size() );
518 for( unsigned ir = 0; ir < verts_n.size(); ir++ )
519 {
520 U( ir, 0 ) = verts_acc2[ir] - verts_acc1[ir];
521 U( ir, 1 ) = verts_acc3[ir] - verts_acc2[ir];
522 res[ir] = -( verts_n[ir] - verts_acc3[ir] );
523 }
524 // U.print();
525 // Vector acc = QR(U).solve(res);
526 Vector acc = solve( U, res );
527 double accsum = acc[0] + acc[1] + 1.0;
528 for( unsigned ir = 0; ir < verts_n.size(); ir++ )
529 {
530 verts_n[ir] = ( verts_acc1[ir] * acc[0] + verts_acc2[ir] * acc[1] + verts_acc3[ir] ) / accsum;
531 }
532
533 memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
534 memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
535 memcpy( &verts_acc3[0], &verts_n[0], nbytes );
536 }
537 else
538 {
539 int offset = 0;
540 for( Range::const_iterator vit = verts.begin(); vit != verts.end(); ++vit, offset += 3 )
541 {
542 // if !fixed
543 if( fix_tag[offset / 3] ) continue;
544
545 CartVect num1 = ( CartVect( &verts_acc3[offset] ) - CartVect( &verts_acc2[offset] ) );
546 CartVect num2 = ( CartVect( &verts_acc3[offset] ) - 2.0 * CartVect( &verts_acc2[offset] ) +
547 CartVect( &verts_acc1[offset] ) );
548
549 num1.scale( num2 );
550 const double mu = num1.length_squared() / num2.length_squared();
551
552 verts_n[offset + 0] =
553 verts_acc3[offset + 0] + mu * ( verts_acc2[offset + 0] - verts_acc3[offset + 0] );
554 verts_n[offset + 1] =
555 verts_acc3[offset + 1] + mu * ( verts_acc2[offset + 1] - verts_acc3[offset + 1] );
556 verts_n[offset + 2] =
557 verts_acc3[offset + 2] + mu * ( verts_acc2[offset + 2] - verts_acc3[offset + 2] );
558 }
559 }
560 }
561 memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
562 memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
563 memcpy( &verts_acc3[0], &verts_n[0], nbytes );
564 }
565
566 // 2b. foreach owned vertex: compute change in coordinate norm
567 for( unsigned v = 0; v < verts.size(); ++v )
568 {
569 double delta = ( CartVect( &verts_n[3 * v] ) - CartVect( &verts_o[3 * v] ) ).length();
570 mxdelta = std::max( delta, mxdelta );
571 EntityHandle vh = verts[v];
572 rval = mb->tag_set_data( errt, &vh, 1, &mxdelta );
573 RC;
574 }
575
576 // global reduce for maximum delta, then report it
577 global_max = mxdelta;
578 #ifdef MOAB_HAVE_MPI
579 if( global_size > 1 ) MPI_Allreduce( &mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, pcomm->comm() );
580 #endif
581
582 if( !( nit % report_its ) )
583 {
584 dbgprint( "\tIterate " << nit << ": Global Max delta = " << global_max << "." );
585 }
586
587 #ifdef WRITE_DEBUG_FILES
588 {
589 // write the tag back onto vertex coordinates
590 rval = mb->set_coords( verts, &verts_n[0] );
591 RC;
592 // output VTK file for debugging purposes
593 std::stringstream sstr;
594 sstr << "LaplacianSmootherIterate_" << nit + 1 << ".h5m";
595 rval = mb->write_file( sstr.str().c_str(), NULL, woptions );
596 RC;
597 }
598 #endif
599
600 #ifdef MOAB_HAVE_MPI
601 // 2c. exchange tags on owned verts
602 if( global_size > 1 )
603 {
604 rval = mb->tag_set_data( vpost, owned_verts, vdata );
605 RC;
606 rval = pcomm->exchange_tags( vpost, shared_owned_verts );
607 RC;
608 }
609 #endif
610
611 if( global_max < rel_eps )
612 break;
613 else
614 {
615 std::copy( verts_n.begin(), verts_n.end(), verts_o.begin() );
616 }
617 }
618
619 // write the tag back onto vertex coordinates
620 rval = mb->set_coords( verts, &verts_n[0] );
621 RC;
622
623 dbgprint( "-- Final iterate error = " << global_max << ".\n" );
624 return MB_SUCCESS;
625 }
626
627 /*
628 Standard Laplacian Smooth Filter
629
630 Additional references: http://www.doc.ic.ac.uk/~gr409/thesis-MSc.pdf
631 */
632 ErrorCode laplacianFilter( Core* mb,
633 moab::ParallelComm* pcomm,
634 moab::Range& verts,
635 int dim,
636 Tag fixed,
637 std::vector< double >& verts_o,
638 std::vector< double >& verts_n,
639 bool use_updated )
640 {
641 ErrorCode rval;
642 std::vector< int > fix_tag( verts.size() );
643 double* data;
644
645 memcpy( &verts_n[0], &verts_o[0], sizeof( double ) * verts_o.size() );
646
647 if( use_updated )
648 data = &verts_n[0];
649 else
650 data = &verts_o[0];
651
652 // filter verts down to owned ones and get fixed tag for them
653 Range owned_verts;
654 if( pcomm->size() > 1 )
655 {
656 #ifdef MOAB_HAVE_MPI
657 rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );
658 if( rval != MB_SUCCESS ) return rval;
659 #endif
660 }
661 else
662 owned_verts = verts;
663
664 rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );
665 RC;
666
667 int vindex = 0;
668 for( Range::const_iterator vit = owned_verts.begin(); vit != owned_verts.end(); ++vit, vindex++ )
669 {
670 // if !fixed
671 if( fix_tag[vindex] ) continue;
672
673 const int index = verts.index( *vit ) * 3;
674
675 moab::Range adjverts, adjelems;
676 // Find the neighboring vertices (1-ring neighborhood)
677 rval = mb->get_adjacencies( &( *vit ), 1, dim, false, adjelems );
678 RC;
679 rval = mb->get_connectivity( adjelems, adjverts );
680 RC;
681 adjverts.erase( *vit );
682
683 const int nadjs = adjverts.size();
684 if( nadjs )
685 {
686 double delta[3] = { 0.0, 0.0, 0.0 };
687
688 // Add the vertices and divide by the number of vertices
689 for( int j = 0; j < nadjs; ++j )
690 {
691 const int joffset = verts.index( adjverts[j] ) * 3;
692 delta[0] += data[joffset + 0];
693 delta[1] += data[joffset + 1];
694 delta[2] += data[joffset + 2];
695 }
696
697 verts_n[index + 0] = delta[0] / nadjs;
698 verts_n[index + 1] = delta[1] / nadjs;
699 verts_n[index + 2] = delta[2] / nadjs;
700 }
701 }
702
703 return MB_SUCCESS;
704 }
705
706 /*
707 HC (Humphrey’s Classes) Smooth Algorithm - Reduces Shrinkage of Laplacian Smoother
708
709 Link:
710 http://informatikbuero.com/downloads/Improved_Laplacian_Smoothing_of_Noisy_Surface_Meshes.pdf
711
712 Where sv - original points
713 pv - previous points,
714 alpha [0..1] influences previous points pv, e.g. 0
715 beta [0..1] e.g. > 0.5
716 */
717 ErrorCode hcFilter( Core* mb,
718 moab::ParallelComm* pcomm,
719 moab::Range& verts,
720 int dim,
721 Tag fixed,
722 std::vector< double >& verts_o,
723 std::vector< double >& verts_n,
724 double alpha,
725 double beta )
726 {
727 ErrorCode rval;
728 std::vector< double > verts_hc( verts_o.size() );
729 std::vector< int > fix_tag( verts.size() );
730
731 // Perform Laplacian Smooth
732 rval = laplacianFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n );
733 RC;
734
735 // Compute Differences
736 for( unsigned index = 0; index < verts_o.size(); ++index )
737 {
738 verts_hc[index] = verts_n[index] - ( alpha * verts_n[index] + ( 1.0 - alpha ) * verts_o[index] );
739 }
740
741 // filter verts down to owned ones and get fixed tag for them
742 Range owned_verts;
743 #ifdef MOAB_HAVE_MPI
744 if( pcomm->size() > 1 )
745 {
746 rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );
747 if( rval != MB_SUCCESS ) return rval;
748 }
749 else
750 #endif
751 owned_verts = verts;
752
753 rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );
754 RC;
755
756 int vindex = 0;
757 for( Range::const_iterator vit = owned_verts.begin(); vit != owned_verts.end(); ++vit, vindex++ )
758 {
759 // if !fixed
760 if( fix_tag[vindex] ) continue;
761
762 const int index = verts.index( *vit ) * 3;
763
764 moab::Range adjverts, adjelems;
765 // Find the neighboring vertices (1-ring neighborhood)
766 rval = mb->get_adjacencies( &( *vit ), 1, dim, false, adjelems );
767 RC;
768 rval = mb->get_connectivity( adjelems, adjverts );
769 RC;
770 adjverts.erase( *vit );
771
772 const int nadjs = adjverts.size();
773 double delta[3] = { 0.0, 0.0, 0.0 };
774
775 for( int j = 0; j < nadjs; ++j )
776 {
777 const int joffset = verts.index( adjverts[j] ) * 3;
778 delta[0] += verts_hc[joffset + 0];
779 delta[1] += verts_hc[joffset + 1];
780 delta[2] += verts_hc[joffset + 2];
781 }
782
783 verts_n[index + 0] -= beta * verts_hc[index + 0] + ( ( 1.0 - beta ) / nadjs ) * delta[0];
784 verts_n[index + 1] -= beta * verts_hc[index + 1] + ( ( 1.0 - beta ) / nadjs ) * delta[1];
785 verts_n[index + 2] -= beta * verts_hc[index + 2] + ( ( 1.0 - beta ) / nadjs ) * delta[2];
786 }
787
788 return MB_SUCCESS;
789 }
1.9.1.