|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
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; }
1.7.6.1 
