Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
LaplacianSmoother.cpp
Go to the documentation of this file.
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 }