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