Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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
#endif
#include "moab/Skinner.hpp"
#include "moab/CN.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 )
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 );
moab::Range& verts,
int dim,
Tag fixed,
std::vector< double >& verts_o,
std::vector< double >& verts_n,
double alpha,
double beta );
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;
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)" ),
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;
}
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
*/
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
*/
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;
}