LaplacianSmoother.cpp

Perform Laplacian relaxation on a mesh and its dual
To run: mpiexec -np <np> LaplacianSmoother [filename]

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

/** @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;
}