LaplacianSmoother.cpp

Go to the documentation of this file.

/** @example LaplacianSmoother.cpp \n
 * \brief Perform Laplacian relaxation on a mesh and its dual \n
 * <b>To run</b>: mpiexec -np <np> LaplacianSmoother [filename]\n
 *
 * Briefly, Laplacian relaxation is a technique to smooth out a mesh. The centroid of each cell is
 * computed from its vertex positions, then vertices are placed at the average of their connected
 * cells' centroids.
 *
 * In the parallel algorithm, an extra ghost layer of cells is exchanged. This allows us to compute
 * the centroids for boundary cells on each processor where they appear; this eliminates the need
 * for one round of data exchange (for those centroids) between processors. New vertex positions
 * must be sent from owning processors to processors sharing those vertices. Convergence is
 * measured as the maximum distance moved by any vertex.
 *
 * In this implementation, a fixed number of iterations is performed. The final mesh is output to
 * 'laplacianfinal.h5m' in the current directory (H5M format must be used since the file is written
 * in parallel).
 *
 * Usage: mpiexec -n 2 valgrind ./LaplacianSmoother -f input/surfrandomtris-64part.h5m -r 2 -p 2 -n
 * 25
 */
 
#include <iostream>
#include <sstream>
#include "moab/Core.hpp"
#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#endif
#include "moab/Skinner.hpp"
#include "moab/CN.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CartVect.hpp"
#include "moab/NestedRefine.hpp"
#include "moab/VerdictWrapper.hpp"
#include "matrix.h"
 
using namespace moab;
using namespace std;
 
#define WRITE_DEBUG_FILES
 
#ifdef MESH_DIR
string test_file_name = string( MESH_DIR ) + string( "/surfrandomtris-4part.h5m" );
#else
string test_file_name = string( "input/surfrandomtris-4part.h5m" );
#endif
 
#define RC MB_CHK_ERR( rval )
#define dbgprint( MSG ) \
 do \
 { \
 if( !global_rank ) std::cerr << ( MSG ) << std::endl; \
 } while( false )
 
ErrorCode perform_laplacian_smoothing( Core* mb,
 Range& cells,
 Range& verts,
 int dim,
 Tag fixed,
 bool use_hc = false,
 bool use_acc = false,
 int acc_method = 1,
 int num_its = 10,
 double rel_eps = 1e-5,
 double alpha = 0.0,
 double beta = 0.5,
 int report_its = 1 );
 
ErrorCode hcFilter( Core* mb,
 moab::ParallelComm* pcomm,
 moab::Range& verts,
 int dim,
 Tag fixed,
 std::vector< double >& verts_o,
 std::vector< double >& verts_n,
 double alpha,
 double beta );
 
ErrorCode laplacianFilter( Core* mb,
 moab::ParallelComm* pcomm,
 moab::Range& verts,
 int dim,
 Tag fixed,
 std::vector< double >& verts_o,
 std::vector< double >& verts_n,
 bool use_updated = true );
 
int main( int argc, char** argv )
{
 int num_its = 10;
 int num_ref = 0;
 int num_dim = 2;
 int report_its = 1;
 int num_degree = 2;
 bool use_hc = false;
 bool use_acc = false;
 int acc_method = 1;
 double alpha = 0.5, beta = 0.0;
 double rel_eps = 1e-5;
 const int nghostrings = 1;
 ProgOptions opts;
 ErrorCode rval;
 std::stringstream sstr;
 int global_rank;
 
 MPI_Init( &argc, &argv );
 MPI_Comm_rank( MPI_COMM_WORLD, &global_rank );
 
 // Decipher program options from user
 opts.addOpt< int >( std::string( "niter,n" ),
 std::string( "Number of Laplacian smoothing iterations (default=10)" ), &num_its );
 opts.addOpt< double >( std::string( "eps,e" ),
 std::string( "Tolerance for the Laplacian smoothing error (default=1e-5)" ), &rel_eps );
 opts.addOpt< double >( std::string( "alpha" ),
 std::string( "Tolerance for the Laplacian smoothing error (default=0.0)" ), &alpha );
 opts.addOpt< double >( std::string( "beta" ),
 std::string( "Tolerance for the Laplacian smoothing error (default=0.5)" ), &beta );
 opts.addOpt< int >( std::string( "dim,d" ), std::string( "Topological dimension of the mesh (default=2)" ),
 &num_dim );
 opts.addOpt< std::string >( std::string( "file,f" ),
 std::string( "Input mesh file to smoothen (default=surfrandomtris-4part.h5m)" ),
 &test_file_name );
 opts.addOpt< int >(
 std::string( "nrefine,r" ),
 std::string( "Number of uniform refinements to perform and apply smoothing cycles (default=1)" ), &num_ref );
 opts.addOpt< int >( std::string( "ndegree,p" ), std::string( "Degree of uniform refinement (default=2)" ),
 &num_degree );
 opts.addOpt< void >( std::string( "humphrey,c" ),
 std::string( "Use Humphrey’s Classes algorithm to reduce shrinkage of "
 "Laplacian smoother (default=false)" ),
 &use_hc );
 opts.addOpt< void >( std::string( "aitken,d" ),
 std::string( "Use Aitken \\delta^2 acceleration to improve convergence of "
 "Lapalace smoothing algorithm (default=false)" ),
 &use_acc );
 opts.addOpt< int >( std::string( "acc,a" ),
 std::string( "Type of vector Aitken process to use for acceleration (default=1)" ),
 &acc_method );
 
 opts.parseCommandLine( argc, argv );
 
 // get MOAB and ParallelComm instances
 Core* mb = new Core;
 if( NULL == mb ) return 1;
 int global_size = 1;
 
 // get the ParallelComm instance
 ParallelComm* pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
 global_size = pcomm->size();
 
 string roptions, woptions;
 if( global_size > 1 )
 { // if reading in parallel, need to tell it how
 sstr.str( "" );
 sstr << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
 "PARALLEL_GHOSTS="
 << num_dim << ".0." << nghostrings << ";DEBUG_IO=0;DEBUG_PIO=0";
 roptions = sstr.str();
 woptions = "PARALLEL=WRITE_PART";
 }
 
 // read the file
 moab::EntityHandle fileset, currset;
 rval = mb->create_meshset( MESHSET_SET, fileset );
 RC;
 currset = fileset;
 rval = mb->load_file( test_file_name.c_str(), &fileset, roptions.c_str() );
 RC;
 
 std::vector< EntityHandle > hsets( num_ref + 1, fileset );
 if( num_ref )
 {
 // Perform uniform refinement of the smoothed mesh
#ifdef MOAB_HAVE_MPI
 NestedRefine* uref = new NestedRefine( mb, pcomm, currset );
#else
 NestedRefine* uref = new NestedRefine( mb, 0, currset );
#endif
 
 std::vector< int > num_degrees( num_ref, num_degree );
 rval = uref->generate_mesh_hierarchy( num_ref, &num_degrees[0], hsets );
 RC;
 
 // Now exchange 1 layer of ghost elements, using vertices as bridge
 // (we could have done this as part of reading process, using the PARALLEL_GHOSTS read
 // option)
 rval = uref->exchange_ghosts( hsets, nghostrings );
 RC;
 
 delete uref;
 }
 
 for( int iref = 0; iref <= num_ref; ++iref )
 {
 
 // specify which set we are currently working on
 currset = hsets[iref];
 
 // make tag to specify fixed vertices, since it's input to the algorithm; use a default
 // value of non-fixed so we only need to set the fixed tag for skin vertices
 Tag fixed;
 int def_val = 0;
 rval = mb->tag_get_handle( "fixed", 1, MB_TYPE_INTEGER, fixed, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );
 RC;
 
 // get all vertices and cells
 Range verts, cells, skin_verts;
 rval = mb->get_entities_by_type( currset, MBVERTEX, verts );
 RC;
 rval = mb->get_entities_by_dimension( currset, num_dim, cells );
 RC;
 dbgprint( "Found " << verts.size() << " vertices and " << cells.size() << " elements" );
 
 // get the skin vertices of those cells and mark them as fixed; we don't want to fix the
 // vertices on a part boundary, but since we exchanged a layer of ghost cells, those
 // vertices aren't on the skin locally ok to mark non-owned skin vertices too, I won't move
 // those anyway use MOAB's skinner class to find the skin
 Skinner skinner( mb );
 rval = skinner.find_skin( currset, cells, true, skin_verts );
 RC; // 'true' param indicates we want vertices back, not cells
 
 std::vector< int > fix_tag( skin_verts.size(), 1 ); // initialized to 1 to indicate fixed
 rval = mb->tag_set_data( fixed, skin_verts, &fix_tag[0] );
 RC;
 // exchange tags on owned verts for fixed points
 if( global_size > 1 )
 {
#ifdef MOAB_HAVE_MPI
 rval = pcomm->exchange_tags( fixed, skin_verts );
 RC;
#endif
 }
 
 // now perform the Laplacian relaxation
 rval = perform_laplacian_smoothing( mb, cells, verts, num_dim, fixed, use_hc, use_acc, acc_method, num_its,
 rel_eps, alpha, beta, report_its );
 RC;
 
 // output file, using parallel write
 sstr.str( "" );
 sstr << "LaplacianSmoother_" << iref << ".h5m";
 rval = mb->write_file( sstr.str().c_str(), "H5M", woptions.c_str(), &currset, 1 );
 RC;
 
 // delete fixed tag, since we created it here
 rval = mb->tag_delete( fixed );
 RC;
 }
 
 delete pcomm;
 // delete MOAB instance
 delete mb;
 
 MPI_Finalize();
 return 0;
}
 
ErrorCode perform_laplacian_smoothing( Core* mb,
 Range& cells,
 Range& verts,
 int dim,
 Tag fixed,
 bool use_hc,
 bool use_acc,
 int acc_method,
 int num_its,
 double rel_eps,
 double alpha,
 double beta,
 int report_its )
{
 ErrorCode rval;
 int global_rank = 0, global_size = 1;
 int nacc = 2; /* nacc_method: 1 = Method 2 from [1], 2 = Method 3 from [1] */
 std::vector< double > verts_acc1, verts_acc2, verts_acc3;
 double rat_theta = rel_eps, rat_alpha = rel_eps, rat_alphaprev = rel_eps;
#ifdef MOAB_HAVE_MPI
 const char* woptions = "PARALLEL=WRITE_PART";
#else
 const char* woptions = "";
#endif
 std::vector< int > fix_tag( verts.size() );
 
 rval = mb->tag_get_data( fixed, verts, &fix_tag[0] );
 RC;
 
#ifdef MOAB_HAVE_MPI
 ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
 global_rank = pcomm->rank();
 global_size = pcomm->size();
#endif
 
 dbgprint( "-- Starting smoothing cycle --" );
 // perform Laplacian relaxation:
 // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
 
 // get all verts coords into tag; don't need to worry about filtering out fixed verts,
 // we'll just be setting to their fixed coords
 std::vector< double > verts_o, verts_n;
 verts_o.resize( 3 * verts.size(), 0.0 );
 verts_n.resize( 3 * verts.size(), 0.0 );
 // std::vector<const void*> vdata(1);
 // vdata[0] = &verts_n[0];
 // const int vdatasize = verts_n.size();
 void* vdata = &verts_n[0];
 rval = mb->get_coords( verts, &verts_o[0] );
 RC;
 const int nbytes = sizeof( double ) * verts_o.size();
 
 Tag errt, vpost;
 double def_val[3] = { 0.0, 0.0, 0.0 };
 rval = mb->tag_get_handle( "error", 1, MB_TYPE_DOUBLE, errt, MB_TAG_CREAT | MB_TAG_DENSE, def_val );
 RC;
 rval = mb->tag_get_handle( "vpos", 3, MB_TYPE_DOUBLE, vpost, MB_TAG_CREAT | MB_TAG_DENSE, def_val );
 RC;
 // rval = mb->tag_set_by_ptr(vpost, verts, vdata); RC;
 
 if( use_acc )
 {
 verts_acc1.resize( verts_o.size(), 0.0 );
 verts_acc2.resize( verts_o.size(), 0.0 );
 verts_acc3.resize( verts_o.size(), 0.0 );
 memcpy( &verts_acc1[0], &verts_o[0], nbytes );
 memcpy( &verts_acc2[0], &verts_o[0], nbytes );
 memcpy( &verts_acc3[0], &verts_o[0], nbytes );
 }
 
 // Filter verts down to owned ones and get fixed tag for them
 Range owned_verts, shared_owned_verts;
 if( global_size > 1 )
 {
#ifdef MOAB_HAVE_MPI
 rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );MB_CHK_ERR( rval );
#endif
 }
 else
 owned_verts = verts;
 
#ifdef MOAB_HAVE_MPI
 // Get shared owned verts, for exchanging tags
 rval = pcomm->get_shared_entities( -1, shared_owned_verts, 0, false, true );MB_CHK_ERR( rval );
 // Workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging
 // tags for all shared entities
 if( shared_owned_verts.empty() ) shared_owned_verts.insert( *verts.begin() );
#endif
 
#ifdef WRITE_DEBUG_FILES
 {
 // output file, using parallel write
 std::stringstream sstr;
 sstr << "LaplacianSmootherIterate_0.h5m";
 rval = mb->write_file( sstr.str().c_str(), NULL, woptions );
 RC;
 }
#endif
 
 bool check_metrics = true;
 VerdictWrapper vw( mb );
 vw.set_size( 1. ); // for relative size measures; maybe it should be user input
 
 double mxdelta = 0.0, global_max = 0.0;
 // 2. for num_its iterations:
 for( int nit = 0; nit < num_its; nit++ )
 {
 
 mxdelta = 0.0;
 
 // 2a. Apply Laplacian Smoothing Filter to Mesh
 if( use_hc )
 {
 rval = hcFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n, alpha, beta );
 RC;
 }
 else
 {
 rval = laplacianFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n );
 RC;
 }
 
 if( check_metrics )
 {
 bool smooth_success = true;
 std::vector< double > coords;
 double jac = 0.0;
 int num_nodes = 0;
 for( unsigned ic = 0; ic < cells.size(); ++ic )
 {
 EntityHandle ecell = cells[ic];
 EntityType etype = mb->type_from_handle( ecell );
 Range everts;
 rval = mb->get_connectivity( &ecell, 1, everts, num_nodes );
 RC;
 coords.resize( num_nodes * 3, 0.0 );
 for( int iv = 0; iv < num_nodes; ++iv )
 {
 const int offset = mb->id_from_handle( everts[iv] ) * 3;
 for( int ii = 0; ii < 3; ++ii )
 coords[iv * 3 + ii] = verts_n[offset + ii];
 }
 rval = vw.quality_measure( ecell, MB_SHAPE, jac, num_nodes, etype, &coords[0] );
 RC;
 if( jac < 1e-8 )
 {
 dbgprint( "Inverted element " << ic << " with jacobian = " << jac << "." );
 smooth_success = false;
 break;
 }
 }
 
 if( !smooth_success )
 {
 verts_o.swap( verts_n );
 dbgprint( "Found element inversions during smoothing. Stopping iterations." );
 break;
 }
 }
 
 if( use_acc )
 {
 
 // if (acc_method < 5 || nit <= 3) {
 // memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
 // memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
 // memcpy( &verts_acc3[0], &verts_n[0], nbytes );
 // }
 
 rat_alphaprev = rat_alpha;
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 rat_alpha =
 std::max( rat_alpha,
 std::abs( ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc2[i] - verts_acc1[i] ) ) /
 ( ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] ) ) );
 }
 rat_theta = std::abs( rat_alpha / rat_alphaprev - 1.0 );
 
 if( nit > 3 && ( nit % nacc ) && rat_theta < 1.0 )
 {
 
 if( acc_method == 1 )
 { /* Method 2 from ACCELERATION OF VECTOR SEQUENCES:
 http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
 double vnorm = 0.0, den, acc_alpha = 0.0, acc_gamma = 0.0;
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 den = ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
 vnorm += den * den;
 acc_alpha += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - verts_acc2[i] );
 acc_gamma += ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] );
 }
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 verts_n[i] = verts_acc2[i] + ( acc_gamma * ( verts_acc3[i] - verts_acc2[i] ) -
 acc_alpha * ( verts_acc2[i] - verts_acc1[i] ) ) /
 vnorm;
 }
 }
 else if( acc_method == 2 )
 { /* Method 3 from ACCELERATION OF VECTOR SEQUENCES:
 http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
 double vnorm = 0.0, num = 0.0, den = 0.0, mu = 0.0;
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 num +=
 ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
 den = ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
 vnorm += den * den;
 }
 mu = num / vnorm;
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 verts_n[i] = verts_acc3[i] + mu * ( verts_acc2[i] - verts_acc3[i] );
 }
 }
 else if( acc_method == 3 )
 { /* Method 5 from ACCELERATION OF VECTOR SEQUENCES:
 http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
 double num = 0.0, den = 0.0, mu = 0.0;
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 num += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc2[i] - verts_acc1[i] );
 den +=
 ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
 }
 mu = num / den;
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 verts_n[i] = verts_acc3[i] - mu * ( verts_acc3[i] - verts_acc2[i] );
 }
 }
 else if( acc_method == 4 )
 { /* Method 8 from ACCELERATION OF VECTOR SEQUENCES:
 http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
 double num = 0.0, den = 0.0, lambda = 0.0;
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 num += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - verts_acc2[i] );
 den += ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] );
 }
 lambda = std::sqrt( num / den );
 for( unsigned i = 0; i < verts_n.size(); ++i )
 {
 verts_n[i] = verts_acc3[i] - lambda / ( lambda - 1.0 ) * ( verts_acc3[i] - verts_acc2[i] );
 }
 }
 else if( acc_method == 5 )
 { /* Minimum polynomial extrapolation of vector sequenes :
 https://en.wikipedia.org/wiki/Minimum_polynomial_extrapolation */
 /* Pseudo-code
 U=x(:,2:end-1)-x(:,1:end-2);
 c=-pinv(U)*(x(:,end)-x(:,end-1));
 c(end+1,1)=1;
 s=(x(:,2:end)*c)/sum(c);
 */
 Matrix U( verts_n.size(), 2 );
 Vector res( verts_n.size() );
 for( unsigned ir = 0; ir < verts_n.size(); ir++ )
 {
 U( ir, 0 ) = verts_acc2[ir] - verts_acc1[ir];
 U( ir, 1 ) = verts_acc3[ir] - verts_acc2[ir];
 res[ir] = -( verts_n[ir] - verts_acc3[ir] );
 }
 // U.print();
 // Vector acc = QR(U).solve(res);
 Vector acc = solve( U, res );
 double accsum = acc[0] + acc[1] + 1.0;
 for( unsigned ir = 0; ir < verts_n.size(); ir++ )
 {
 verts_n[ir] = ( verts_acc1[ir] * acc[0] + verts_acc2[ir] * acc[1] + verts_acc3[ir] ) / accsum;
 }
 
 memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
 memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
 memcpy( &verts_acc3[0], &verts_n[0], nbytes );
 }
 else
 {
 int offset = 0;
 for( Range::const_iterator vit = verts.begin(); vit != verts.end(); ++vit, offset += 3 )
 {
 // if !fixed
 if( fix_tag[offset / 3] ) continue;
 
 CartVect num1 = ( CartVect( &verts_acc3[offset] ) - CartVect( &verts_acc2[offset] ) );
 CartVect num2 = ( CartVect( &verts_acc3[offset] ) - 2.0 * CartVect( &verts_acc2[offset] ) +
 CartVect( &verts_acc1[offset] ) );
 
 num1.scale( num2 );
 const double mu = num1.length_squared() / num2.length_squared();
 
 verts_n[offset + 0] =
 verts_acc3[offset + 0] + mu * ( verts_acc2[offset + 0] - verts_acc3[offset + 0] );
 verts_n[offset + 1] =
 verts_acc3[offset + 1] + mu * ( verts_acc2[offset + 1] - verts_acc3[offset + 1] );
 verts_n[offset + 2] =
 verts_acc3[offset + 2] + mu * ( verts_acc2[offset + 2] - verts_acc3[offset + 2] );
 }
 }
 }
 memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
 memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
 memcpy( &verts_acc3[0], &verts_n[0], nbytes );
 }
 
 // 2b. foreach owned vertex: compute change in coordinate norm
 for( unsigned v = 0; v < verts.size(); ++v )
 {
 double delta = ( CartVect( &verts_n[3 * v] ) - CartVect( &verts_o[3 * v] ) ).length();
 mxdelta = std::max( delta, mxdelta );
 EntityHandle vh = verts[v];
 rval = mb->tag_set_data( errt, &vh, 1, &mxdelta );
 RC;
 }
 
 // global reduce for maximum delta, then report it
 global_max = mxdelta;
#ifdef MOAB_HAVE_MPI
 if( global_size > 1 ) MPI_Allreduce( &mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, pcomm->comm() );
#endif
 
 if( !( nit % report_its ) )
 {
 dbgprint( "\tIterate " << nit << ": Global Max delta = " << global_max << "." );
 }
 
#ifdef WRITE_DEBUG_FILES
 {
 // write the tag back onto vertex coordinates
 rval = mb->set_coords( verts, &verts_n[0] );
 RC;
 // output VTK file for debugging purposes
 std::stringstream sstr;
 sstr << "LaplacianSmootherIterate_" << nit + 1 << ".h5m";
 rval = mb->write_file( sstr.str().c_str(), NULL, woptions );
 RC;
 }
#endif
 
#ifdef MOAB_HAVE_MPI
 // 2c. exchange tags on owned verts
 if( global_size > 1 )
 {
 rval = mb->tag_set_data( vpost, owned_verts, vdata );
 RC;
 rval = pcomm->exchange_tags( vpost, shared_owned_verts );
 RC;
 }
#endif
 
 if( global_max < rel_eps )
 break;
 else
 {
 std::copy( verts_n.begin(), verts_n.end(), verts_o.begin() );
 }
 }
 
 // write the tag back onto vertex coordinates
 rval = mb->set_coords( verts, &verts_n[0] );
 RC;
 
 dbgprint( "-- Final iterate error = " << global_max << ".\n" );
 return MB_SUCCESS;
}
 
/*
 Standard Laplacian Smooth Filter
 
 Additional references: http://www.doc.ic.ac.uk/~gr409/thesis-MSc.pdf
*/
ErrorCode laplacianFilter( Core* mb,
 moab::ParallelComm* pcomm,
 moab::Range& verts,
 int dim,
 Tag fixed,
 std::vector< double >& verts_o,
 std::vector< double >& verts_n,
 bool use_updated )
{
 ErrorCode rval;
 std::vector< int > fix_tag( verts.size() );
 double* data;
 
 memcpy( &verts_n[0], &verts_o[0], sizeof( double ) * verts_o.size() );
 
 if( use_updated )
 data = &verts_n[0];
 else
 data = &verts_o[0];
 
 // filter verts down to owned ones and get fixed tag for them
 Range owned_verts;
 if( pcomm->size() > 1 )
 {
#ifdef MOAB_HAVE_MPI
 rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );
 if( rval != MB_SUCCESS ) return rval;
#endif
 }
 else
 owned_verts = verts;
 
 rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );
 RC;
 
 int vindex = 0;
 for( Range::const_iterator vit = owned_verts.begin(); vit != owned_verts.end(); ++vit, vindex++ )
 {
 // if !fixed
 if( fix_tag[vindex] ) continue;
 
 const int index = verts.index( *vit ) * 3;
 
 moab::Range adjverts, adjelems;
 // Find the neighboring vertices (1-ring neighborhood)
 rval = mb->get_adjacencies( &( *vit ), 1, dim, false, adjelems );
 RC;
 rval = mb->get_connectivity( adjelems, adjverts );
 RC;
 adjverts.erase( *vit );
 
 const int nadjs = adjverts.size();
 if( nadjs )
 {
 double delta[3] = { 0.0, 0.0, 0.0 };
 
 // Add the vertices and divide by the number of vertices
 for( int j = 0; j < nadjs; ++j )
 {
 const int joffset = verts.index( adjverts[j] ) * 3;
 delta[0] += data[joffset + 0];
 delta[1] += data[joffset + 1];
 delta[2] += data[joffset + 2];
 }
 
 verts_n[index + 0] = delta[0] / nadjs;
 verts_n[index + 1] = delta[1] / nadjs;
 verts_n[index + 2] = delta[2] / nadjs;
 }
 }
 
 return MB_SUCCESS;
}
 
/*
 HC (Humphrey’s Classes) Smooth Algorithm - Reduces Shrinkage of Laplacian Smoother
 
 Link:
 http://informatikbuero.com/downloads/Improved_Laplacian_Smoothing_of_Noisy_Surface_Meshes.pdf
 
 Where sv - original points
 pv - previous points,
 alpha [0..1] influences previous points pv, e.g. 0
 beta [0..1] e.g. > 0.5
*/
ErrorCode hcFilter( Core* mb,
 moab::ParallelComm* pcomm,
 moab::Range& verts,
 int dim,
 Tag fixed,
 std::vector< double >& verts_o,
 std::vector< double >& verts_n,
 double alpha,
 double beta )
{
 ErrorCode rval;
 std::vector< double > verts_hc( verts_o.size() );
 std::vector< int > fix_tag( verts.size() );
 
 // Perform Laplacian Smooth
 rval = laplacianFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n );
 RC;
 
 // Compute Differences
 for( unsigned index = 0; index < verts_o.size(); ++index )
 {
 verts_hc[index] = verts_n[index] - ( alpha * verts_n[index] + ( 1.0 - alpha ) * verts_o[index] );
 }
 
 // filter verts down to owned ones and get fixed tag for them
 Range owned_verts;
#ifdef MOAB_HAVE_MPI
 if( pcomm->size() > 1 )
 {
 rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );
 if( rval != MB_SUCCESS ) return rval;
 }
 else
#endif
 owned_verts = verts;
 
 rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );
 RC;
 
 int vindex = 0;
 for( Range::const_iterator vit = owned_verts.begin(); vit != owned_verts.end(); ++vit, vindex++ )
 {
 // if !fixed
 if( fix_tag[vindex] ) continue;
 
 const int index = verts.index( *vit ) * 3;
 
 moab::Range adjverts, adjelems;
 // Find the neighboring vertices (1-ring neighborhood)
 rval = mb->get_adjacencies( &( *vit ), 1, dim, false, adjelems );
 RC;
 rval = mb->get_connectivity( adjelems, adjverts );
 RC;
 adjverts.erase( *vit );
 
 const int nadjs = adjverts.size();
 double delta[3] = { 0.0, 0.0, 0.0 };
 
 for( int j = 0; j < nadjs; ++j )
 {
 const int joffset = verts.index( adjverts[j] ) * 3;
 delta[0] += verts_hc[joffset + 0];
 delta[1] += verts_hc[joffset + 1];
 delta[2] += verts_hc[joffset + 2];
 }
 
 verts_n[index + 0] -= beta * verts_hc[index + 0] + ( ( 1.0 - beta ) / nadjs ) * delta[0];
 verts_n[index + 1] -= beta * verts_hc[index + 1] + ( ( 1.0 - beta ) / nadjs ) * delta[1];
 verts_n[index + 2] -= beta * verts_hc[index + 2] + ( ( 1.0 - beta ) / nadjs ) * delta[2];
 }
 
 return MB_SUCCESS;
}