MOAB: Mesh Oriented datABase  (version 5.2.1)
LaplacianSmoother.cpp
Go to the documentation of this file.
00001 /** @example LaplacianSmoother.cpp \n
00002  * \brief Perform Laplacian relaxation on a mesh and its dual \n
00003  * <b>To run</b>: mpiexec -np <np> LaplacianSmoother [filename]\n
00004  *
00005  * Briefly, Laplacian relaxation is a technique to smooth out a mesh.  The centroid of each cell is
00006  * computed from its vertex positions, then vertices are placed at the average of their connected
00007  * cells' centroids.
00008  *
00009  * In the parallel algorithm, an extra ghost layer of cells is exchanged.  This allows us to compute
00010  * the centroids for boundary cells on each processor where they appear; this eliminates the need
00011  * for one round of data exchange (for those centroids) between processors.  New vertex positions
00012  * must be sent from owning processors to processors sharing those vertices.  Convergence is
00013  * measured as the maximum distance moved by any vertex.
00014  *
00015  * In this implementation, a fixed number of iterations is performed.  The final mesh is output to
00016  * 'laplacianfinal.h5m' in the current directory (H5M format must be used since the file is written
00017  * in parallel).
00018  *
00019  * Usage: mpiexec -n 2 valgrind ./LaplacianSmoother -f input/surfrandomtris-64part.h5m -r 2 -p 2 -n
00020  * 25
00021  */
00022 
00023 #include <iostream>
00024 #include <sstream>
00025 #include "moab/Core.hpp"
00026 #ifdef MOAB_HAVE_MPI
00027 #include "moab/ParallelComm.hpp"
00028 #include "MBParallelConventions.h"
00029 #endif
00030 #include "moab/Skinner.hpp"
00031 #include "moab/CN.hpp"
00032 #include "moab/ProgOptions.hpp"
00033 #include "moab/CartVect.hpp"
00034 #include "moab/NestedRefine.hpp"
00035 #include "moab/VerdictWrapper.hpp"
00036 #include "matrix.h"
00037 
00038 using namespace moab;
00039 using namespace std;
00040 
00041 #define WRITE_DEBUG_FILES
00042 
00043 #ifdef MESH_DIR
00044 string test_file_name = string( MESH_DIR ) + string( "/surfrandomtris-4part.h5m" );
00045 #else
00046 string test_file_name = string( "input/surfrandomtris-4part.h5m" );
00047 #endif
00048 
00049 #define RC MB_CHK_ERR( rval )
00050 #define dbgprint( MSG )                                   \
00051     do                                                    \
00052     {                                                     \
00053         if( !global_rank ) std::cerr << MSG << std::endl; \
00054     } while( false )
00055 
00056 ErrorCode perform_laplacian_smoothing( Core* mb, Range& cells, Range& verts, int dim, Tag fixed, bool use_hc = false,
00057                                        bool use_acc = false, int acc_method = 1, int num_its = 10,
00058                                        double rel_eps = 1e-5, double alpha = 0.0, double beta = 0.5,
00059                                        int report_its = 1 );
00060 
00061 ErrorCode hcFilter( Core* mb, moab::ParallelComm* pcomm, moab::Range& verts, int dim, Tag fixed,
00062                     std::vector< double >& verts_o, std::vector< double >& verts_n, double alpha, double beta );
00063 
00064 ErrorCode laplacianFilter( Core* mb, moab::ParallelComm* pcomm, moab::Range& verts, int dim, Tag fixed,
00065                            std::vector< double >& verts_o, std::vector< double >& verts_n, bool use_updated = true );
00066 
00067 int main( int argc, char** argv )
00068 {
00069     int num_its    = 10;
00070     int num_ref    = 0;
00071     int num_dim    = 2;
00072     int report_its = 1;
00073     int num_degree = 2;
00074     bool use_hc    = false;
00075     bool use_acc   = false;
00076     int acc_method = 1;
00077     double alpha = 0.5, beta = 0.0;
00078     double rel_eps        = 1e-5;
00079     const int nghostrings = 1;
00080     ProgOptions opts;
00081     ErrorCode rval;
00082     std::stringstream sstr;
00083     int global_rank;
00084 
00085     MPI_Init( &argc, &argv );
00086     MPI_Comm_rank( MPI_COMM_WORLD, &global_rank );
00087 
00088     // Decipher program options from user
00089     opts.addOpt< int >( std::string( "niter,n" ),
00090                         std::string( "Number of Laplacian smoothing iterations (default=10)" ), &num_its );
00091     opts.addOpt< double >( std::string( "eps,e" ),
00092                            std::string( "Tolerance for the Laplacian smoothing error (default=1e-5)" ), &rel_eps );
00093     opts.addOpt< double >( std::string( "alpha" ),
00094                            std::string( "Tolerance for the Laplacian smoothing error (default=0.0)" ), &alpha );
00095     opts.addOpt< double >( std::string( "beta" ),
00096                            std::string( "Tolerance for the Laplacian smoothing error (default=0.5)" ), &beta );
00097     opts.addOpt< int >( std::string( "dim,d" ), std::string( "Topological dimension of the mesh (default=2)" ),
00098                         &num_dim );
00099     opts.addOpt< std::string >( std::string( "file,f" ),
00100                                 std::string( "Input mesh file to smoothen (default=surfrandomtris-4part.h5m)" ),
00101                                 &test_file_name );
00102     opts.addOpt< int >(
00103         std::string( "nrefine,r" ),
00104         std::string( "Number of uniform refinements to perform and apply smoothing cycles (default=1)" ), &num_ref );
00105     opts.addOpt< int >( std::string( "ndegree,p" ), std::string( "Degree of uniform refinement (default=2)" ),
00106                         &num_degree );
00107     opts.addOpt< void >( std::string( "humphrey,c" ),
00108                          std::string( "Use Humphrey’s Classes algorithm to reduce shrinkage of "
00109                                       "Laplacian smoother (default=false)" ),
00110                          &use_hc );
00111     opts.addOpt< void >( std::string( "aitken,d" ),
00112                          std::string( "Use Aitken \\delta^2 acceleration to improve convergence of "
00113                                       "Lapalace smoothing algorithm (default=false)" ),
00114                          &use_acc );
00115     opts.addOpt< int >( std::string( "acc,a" ),
00116                         std::string( "Type of vector Aitken process to use for acceleration (default=1)" ),
00117                         &acc_method );
00118 
00119     opts.parseCommandLine( argc, argv );
00120 
00121     // get MOAB and ParallelComm instances
00122     Core* mb = new Core;
00123     if( NULL == mb ) return 1;
00124     int global_size = 1;
00125 
00126     // get the ParallelComm instance
00127     ParallelComm* pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
00128     global_size         = pcomm->size();
00129 
00130     string roptions, woptions;
00131     if( global_size > 1 )
00132     {  // if reading in parallel, need to tell it how
00133         sstr.str( "" );
00134         sstr << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
00135                 "PARALLEL_GHOSTS="
00136              << num_dim << ".0." << nghostrings << ";DEBUG_IO=0;DEBUG_PIO=0";
00137         roptions = sstr.str();
00138         woptions = "PARALLEL=WRITE_PART";
00139     }
00140 
00141     // read the file
00142     moab::EntityHandle fileset, currset;
00143     rval = mb->create_meshset( MESHSET_SET, fileset );
00144     RC;
00145     currset = fileset;
00146     rval    = mb->load_file( test_file_name.c_str(), &fileset, roptions.c_str() );
00147     RC;
00148 
00149     std::vector< EntityHandle > hsets( num_ref + 1, fileset );
00150     if( num_ref )
00151     {
00152         // Perform uniform refinement of the smoothed mesh
00153 #ifdef MOAB_HAVE_MPI
00154         NestedRefine* uref = new NestedRefine( mb, pcomm, currset );
00155 #else
00156         NestedRefine* uref = new NestedRefine( mb, 0, currset );
00157 #endif
00158 
00159         std::vector< int > num_degrees( num_ref, num_degree );
00160         rval = uref->generate_mesh_hierarchy( num_ref, &num_degrees[0], hsets );
00161         RC;
00162 
00163         // Now exchange 1 layer of ghost elements, using vertices as bridge
00164         // (we could have done this as part of reading process, using the PARALLEL_GHOSTS read
00165         // option)
00166         rval = uref->exchange_ghosts( hsets, nghostrings );
00167         RC;
00168 
00169         delete uref;
00170     }
00171 
00172     for( int iref = 0; iref <= num_ref; ++iref )
00173     {
00174 
00175         // specify which set we are currently working on
00176         currset = hsets[iref];
00177 
00178         // make tag to specify fixed vertices, since it's input to the algorithm; use a default
00179         // value of non-fixed so we only need to set the fixed tag for skin vertices
00180         Tag fixed;
00181         int def_val = 0;
00182         rval        = mb->tag_get_handle( "fixed", 1, MB_TYPE_INTEGER, fixed, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );
00183         RC;
00184 
00185         // get all vertices and cells
00186         Range verts, cells, skin_verts;
00187         rval = mb->get_entities_by_type( currset, MBVERTEX, verts );
00188         RC;
00189         rval = mb->get_entities_by_dimension( currset, num_dim, cells );
00190         RC;
00191         dbgprint( "Found " << verts.size() << " vertices and " << cells.size() << " elements" );
00192 
00193         // get the skin vertices of those cells and mark them as fixed; we don't want to fix the
00194         // vertices on a part boundary, but since we exchanged a layer of ghost cells, those
00195         // vertices aren't on the skin locally ok to mark non-owned skin vertices too, I won't move
00196         // those anyway use MOAB's skinner class to find the skin
00197         Skinner skinner( mb );
00198         rval = skinner.find_skin( currset, cells, true, skin_verts );
00199         RC;  // 'true' param indicates we want vertices back, not cells
00200 
00201         std::vector< int > fix_tag( skin_verts.size(), 1 );  // initialized to 1 to indicate fixed
00202         rval = mb->tag_set_data( fixed, skin_verts, &fix_tag[0] );
00203         RC;
00204         // exchange tags on owned verts for fixed points
00205         if( global_size > 1 )
00206         {
00207 #ifdef MOAB_HAVE_MPI
00208             rval = pcomm->exchange_tags( fixed, skin_verts );
00209             RC;
00210 #endif
00211         }
00212 
00213         // now perform the Laplacian relaxation
00214         rval = perform_laplacian_smoothing( mb, cells, verts, num_dim, fixed, use_hc, use_acc, acc_method, num_its,
00215                                             rel_eps, alpha, beta, report_its );
00216         RC;
00217 
00218         // output file, using parallel write
00219         sstr.str( "" );
00220         sstr << "LaplacianSmoother_" << iref << ".h5m";
00221         rval = mb->write_file( sstr.str().c_str(), "H5M", woptions.c_str(), &currset, 1 );
00222         RC;
00223 
00224         // delete fixed tag, since we created it here
00225         rval = mb->tag_delete( fixed );
00226         RC;
00227     }
00228 
00229     delete pcomm;
00230     // delete MOAB instance
00231     delete mb;
00232 
00233     MPI_Finalize();
00234     return 0;
00235 }
00236 
00237 ErrorCode perform_laplacian_smoothing( Core* mb, Range& cells, Range& verts, int dim, Tag fixed, bool use_hc,
00238                                        bool use_acc, int acc_method, int num_its, double rel_eps, double alpha,
00239                                        double beta, int report_its )
00240 {
00241     ErrorCode rval;
00242     int global_rank = 0, global_size = 1;
00243     int nacc = 2; /* nacc_method: 1 = Method 2 from [1], 2 = Method 3 from [1] */
00244     std::vector< double > verts_acc1, verts_acc2, verts_acc3;
00245     double rat_theta = rel_eps, rat_alpha = rel_eps, rat_alphaprev = rel_eps;
00246 #ifdef MOAB_HAVE_MPI
00247     const char* woptions = "PARALLEL=WRITE_PART";
00248 #else
00249     const char* woptions = "";
00250 #endif
00251     std::vector< int > fix_tag( verts.size() );
00252 
00253     rval = mb->tag_get_data( fixed, verts, &fix_tag[0] );
00254     RC;
00255 
00256 #ifdef MOAB_HAVE_MPI
00257     ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
00258     global_rank         = pcomm->rank();
00259     global_size         = pcomm->size();
00260 #endif
00261 
00262     dbgprint( "-- Starting smoothing cycle --" );
00263     // perform Laplacian relaxation:
00264     // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
00265 
00266     // get all verts coords into tag; don't need to worry about filtering out fixed verts,
00267     // we'll just be setting to their fixed coords
00268     std::vector< double > verts_o, verts_n;
00269     verts_o.resize( 3 * verts.size(), 0.0 );
00270     verts_n.resize( 3 * verts.size(), 0.0 );
00271     // std::vector<const void*> vdata(1);
00272     // vdata[0] = &verts_n[0];
00273     // const int vdatasize = verts_n.size();
00274     void* vdata = &verts_n[0];
00275     rval        = mb->get_coords( verts, &verts_o[0] );
00276     RC;
00277     const int nbytes = sizeof( double ) * verts_o.size();
00278 
00279     Tag errt, vpost;
00280     double def_val[3] = { 0.0, 0.0, 0.0 };
00281     rval              = mb->tag_get_handle( "error", 1, MB_TYPE_DOUBLE, errt, MB_TAG_CREAT | MB_TAG_DENSE, def_val );
00282     RC;
00283     rval = mb->tag_get_handle( "vpos", 3, MB_TYPE_DOUBLE, vpost, MB_TAG_CREAT | MB_TAG_DENSE, def_val );
00284     RC;
00285     // rval = mb->tag_set_by_ptr(vpost, verts, vdata); RC;
00286 
00287     if( use_acc )
00288     {
00289         verts_acc1.resize( verts_o.size(), 0.0 );
00290         verts_acc2.resize( verts_o.size(), 0.0 );
00291         verts_acc3.resize( verts_o.size(), 0.0 );
00292         memcpy( &verts_acc1[0], &verts_o[0], nbytes );
00293         memcpy( &verts_acc2[0], &verts_o[0], nbytes );
00294         memcpy( &verts_acc3[0], &verts_o[0], nbytes );
00295     }
00296 
00297     // Filter verts down to owned ones and get fixed tag for them
00298     Range owned_verts, shared_owned_verts;
00299     if( global_size > 1 )
00300     {
00301 #ifdef MOAB_HAVE_MPI
00302         rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );MB_CHK_ERR( rval );
00303 #endif
00304     }
00305     else
00306         owned_verts = verts;
00307 
00308 #ifdef MOAB_HAVE_MPI
00309     // Get shared owned verts, for exchanging tags
00310     rval = pcomm->get_shared_entities( -1, shared_owned_verts, 0, false, true );MB_CHK_ERR( rval );
00311     // Workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging
00312     // tags for all shared entities
00313     if( shared_owned_verts.empty() ) shared_owned_verts.insert( *verts.begin() );
00314 #endif
00315 
00316 #ifdef WRITE_DEBUG_FILES
00317     {
00318         // output file, using parallel write
00319         std::stringstream sstr;
00320         sstr << "LaplacianSmootherIterate_0.h5m";
00321         rval = mb->write_file( sstr.str().c_str(), NULL, woptions );
00322         RC;
00323     }
00324 #endif
00325 
00326     bool check_metrics = true;
00327     VerdictWrapper vw( mb );
00328     vw.set_size( 1. );  // for relative size measures; maybe it should be user input
00329 
00330     double mxdelta = 0.0, global_max = 0.0;
00331     // 2. for num_its iterations:
00332     for( int nit = 0; nit < num_its; nit++ )
00333     {
00334 
00335         mxdelta = 0.0;
00336 
00337         // 2a. Apply Laplacian Smoothing Filter to Mesh
00338         if( use_hc )
00339         {
00340             rval = hcFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n, alpha, beta );
00341             RC;
00342         }
00343         else
00344         {
00345             rval = laplacianFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n );
00346             RC;
00347         }
00348 
00349         if( check_metrics )
00350         {
00351             bool smooth_success = true;
00352             std::vector< double > coords;
00353             double jac    = 0.0;
00354             int num_nodes = 0;
00355             for( unsigned ic = 0; ic < cells.size(); ++ic )
00356             {
00357                 EntityHandle ecell = cells[ic];
00358                 EntityType etype   = mb->type_from_handle( ecell );
00359                 Range everts;
00360                 rval = mb->get_connectivity( &ecell, 1, everts, num_nodes );
00361                 RC;
00362                 coords.resize( num_nodes * 3, 0.0 );
00363                 for( int iv = 0; iv < num_nodes; ++iv )
00364                 {
00365                     const int offset = mb->id_from_handle( everts[iv] ) * 3;
00366                     for( int ii = 0; ii < 3; ++ii )
00367                         coords[iv * 3 + ii] = verts_n[offset + ii];
00368                 }
00369                 rval = vw.quality_measure( ecell, MB_SHAPE, jac, num_nodes, etype, &coords[0] );
00370                 RC;
00371                 if( jac < 1e-8 )
00372                 {
00373                     dbgprint( "Inverted element " << ic << " with jacobian = " << jac << "." );
00374                     smooth_success = false;
00375                     break;
00376                 }
00377             }
00378 
00379             if( !smooth_success )
00380             {
00381                 verts_o.swap( verts_n );
00382                 dbgprint( "Found element inversions during smoothing. Stopping iterations." );
00383                 break;
00384             }
00385         }
00386 
00387         if( use_acc )
00388         {
00389 
00390             // if (acc_method < 5 || nit <= 3) {
00391             //   memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
00392             //   memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
00393             //   memcpy( &verts_acc3[0], &verts_n[0],    nbytes );
00394             // }
00395 
00396             rat_alphaprev = rat_alpha;
00397             for( unsigned i = 0; i < verts_n.size(); ++i )
00398             {
00399                 rat_alpha =
00400                     std::max( rat_alpha,
00401                               std::abs( ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc2[i] - verts_acc1[i] ) ) /
00402                                   ( ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] ) ) );
00403             }
00404             rat_theta = std::abs( rat_alpha / rat_alphaprev - 1.0 );
00405 
00406             if( nit > 3 && ( nit % nacc ) && rat_theta < 1.0 )
00407             {
00408 
00409                 if( acc_method == 1 )
00410                 { /* Method 2 from ACCELERATION OF VECTOR SEQUENCES:
00411                      http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
00412                     double vnorm = 0.0, den, acc_alpha = 0.0, acc_gamma = 0.0;
00413                     for( unsigned i = 0; i < verts_n.size(); ++i )
00414                     {
00415                         den = ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
00416                         vnorm += den * den;
00417                         acc_alpha += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - verts_acc2[i] );
00418                         acc_gamma += ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] );
00419                     }
00420                     for( unsigned i = 0; i < verts_n.size(); ++i )
00421                     {
00422                         verts_n[i] = verts_acc2[i] + ( acc_gamma * ( verts_acc3[i] - verts_acc2[i] ) -
00423                                                        acc_alpha * ( verts_acc2[i] - verts_acc1[i] ) ) /
00424                                                          vnorm;
00425                     }
00426                 }
00427                 else if( acc_method == 2 )
00428                 { /* Method 3 from ACCELERATION OF VECTOR SEQUENCES:
00429                      http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
00430                     double vnorm = 0.0, num = 0.0, den = 0.0, mu = 0.0;
00431                     for( unsigned i = 0; i < verts_n.size(); ++i )
00432                     {
00433                         num +=
00434                             ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
00435                         den = ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
00436                         vnorm += den * den;
00437                     }
00438                     mu = num / vnorm;
00439                     for( unsigned i = 0; i < verts_n.size(); ++i )
00440                     {
00441                         verts_n[i] = verts_acc3[i] + mu * ( verts_acc2[i] - verts_acc3[i] );
00442                     }
00443                 }
00444                 else if( acc_method == 3 )
00445                 { /* Method 5 from ACCELERATION OF VECTOR SEQUENCES:
00446                      http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
00447                     double num = 0.0, den = 0.0, mu = 0.0;
00448                     for( unsigned i = 0; i < verts_n.size(); ++i )
00449                     {
00450                         num += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc2[i] - verts_acc1[i] );
00451                         den +=
00452                             ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc3[i] - 2.0 * verts_acc2[i] + verts_acc1[i] );
00453                     }
00454                     mu = num / den;
00455                     for( unsigned i = 0; i < verts_n.size(); ++i )
00456                     {
00457                         verts_n[i] = verts_acc3[i] - mu * ( verts_acc3[i] - verts_acc2[i] );
00458                     }
00459                 }
00460                 else if( acc_method == 4 )
00461                 { /* Method 8 from ACCELERATION OF VECTOR SEQUENCES:
00462                      http://onlinelibrary.wiley.com/doi/10.1002/cnm.1630020409/pdf */
00463                     double num = 0.0, den = 0.0, lambda = 0.0;
00464                     for( unsigned i = 0; i < verts_n.size(); ++i )
00465                     {
00466                         num += ( verts_acc3[i] - verts_acc2[i] ) * ( verts_acc3[i] - verts_acc2[i] );
00467                         den += ( verts_acc2[i] - verts_acc1[i] ) * ( verts_acc2[i] - verts_acc1[i] );
00468                     }
00469                     lambda = std::sqrt( num / den );
00470                     for( unsigned i = 0; i < verts_n.size(); ++i )
00471                     {
00472                         verts_n[i] = verts_acc3[i] - lambda / ( lambda - 1.0 ) * ( verts_acc3[i] - verts_acc2[i] );
00473                     }
00474                 }
00475                 else if( acc_method == 5 )
00476                 { /* Minimum polynomial extrapolation of vector sequenes :
00477                      https://en.wikipedia.org/wiki/Minimum_polynomial_extrapolation */
00478                     /* Pseudo-code
00479                           U=x(:,2:end-1)-x(:,1:end-2);
00480                           c=-pinv(U)*(x(:,end)-x(:,end-1));
00481                           c(end+1,1)=1;
00482                           s=(x(:,2:end)*c)/sum(c);
00483                     */
00484                     Matrix U( verts_n.size(), 2 );
00485                     Vector res( verts_n.size() );
00486                     for( unsigned ir = 0; ir < verts_n.size(); ir++ )
00487                     {
00488                         U( ir, 0 ) = verts_acc2[ir] - verts_acc1[ir];
00489                         U( ir, 1 ) = verts_acc3[ir] - verts_acc2[ir];
00490                         res[ir]    = -( verts_n[ir] - verts_acc3[ir] );
00491                     }
00492                     // U.print();
00493                     // Vector acc = QR(U).solve(res);
00494                     Vector acc    = solve( U, res );
00495                     double accsum = acc[0] + acc[1] + 1.0;
00496                     for( unsigned ir = 0; ir < verts_n.size(); ir++ )
00497                     {
00498                         verts_n[ir] = ( verts_acc1[ir] * acc[0] + verts_acc2[ir] * acc[1] + verts_acc3[ir] ) / accsum;
00499                     }
00500 
00501                     memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
00502                     memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
00503                     memcpy( &verts_acc3[0], &verts_n[0], nbytes );
00504                 }
00505                 else
00506                 {
00507                     int offset = 0;
00508                     for( Range::const_iterator vit = verts.begin(); vit != verts.end(); ++vit, offset += 3 )
00509                     {
00510                         // if !fixed
00511                         if( fix_tag[offset / 3] ) continue;
00512 
00513                         CartVect num1 = ( CartVect( &verts_acc3[offset] ) - CartVect( &verts_acc2[offset] ) );
00514                         CartVect num2 = ( CartVect( &verts_acc3[offset] ) - 2.0 * CartVect( &verts_acc2[offset] ) +
00515                                           CartVect( &verts_acc1[offset] ) );
00516 
00517                         num1.scale( num2 );
00518                         const double mu = num1.length_squared() / num2.length_squared();
00519 
00520                         verts_n[offset + 0] =
00521                             verts_acc3[offset + 0] + mu * ( verts_acc2[offset + 0] - verts_acc3[offset + 0] );
00522                         verts_n[offset + 1] =
00523                             verts_acc3[offset + 1] + mu * ( verts_acc2[offset + 1] - verts_acc3[offset + 1] );
00524                         verts_n[offset + 2] =
00525                             verts_acc3[offset + 2] + mu * ( verts_acc2[offset + 2] - verts_acc3[offset + 2] );
00526                     }
00527                 }
00528             }
00529             memcpy( &verts_acc1[0], &verts_acc2[0], nbytes );
00530             memcpy( &verts_acc2[0], &verts_acc3[0], nbytes );
00531             memcpy( &verts_acc3[0], &verts_n[0], nbytes );
00532         }
00533 
00534         // 2b. foreach owned vertex: compute change in coordinate norm
00535         for( unsigned v = 0; v < verts.size(); ++v )
00536         {
00537             double delta    = ( CartVect( &verts_n[3 * v] ) - CartVect( &verts_o[3 * v] ) ).length();
00538             mxdelta         = std::max( delta, mxdelta );
00539             EntityHandle vh = verts[v];
00540             rval            = mb->tag_set_data( errt, &vh, 1, &mxdelta );
00541             RC;
00542         }
00543 
00544         // global reduce for maximum delta, then report it
00545         global_max = mxdelta;
00546 #ifdef MOAB_HAVE_MPI
00547         if( global_size > 1 ) MPI_Allreduce( &mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, pcomm->comm() );
00548 #endif
00549 
00550         if( !( nit % report_its ) ) { dbgprint( "\tIterate " << nit << ": Global Max delta = " << global_max << "." ); }
00551 
00552 #ifdef WRITE_DEBUG_FILES
00553         {
00554             // write the tag back onto vertex coordinates
00555             rval = mb->set_coords( verts, &verts_n[0] );
00556             RC;
00557             // output VTK file for debugging purposes
00558             std::stringstream sstr;
00559             sstr << "LaplacianSmootherIterate_" << nit + 1 << ".h5m";
00560             rval = mb->write_file( sstr.str().c_str(), NULL, woptions );
00561             RC;
00562         }
00563 #endif
00564 
00565 #ifdef MOAB_HAVE_MPI
00566         // 2c. exchange tags on owned verts
00567         if( global_size > 1 )
00568         {
00569             rval = mb->tag_set_data( vpost, owned_verts, vdata );
00570             RC;
00571             rval = pcomm->exchange_tags( vpost, shared_owned_verts );
00572             RC;
00573         }
00574 #endif
00575 
00576         if( global_max < rel_eps )
00577             break;
00578         else
00579         {
00580             std::copy( verts_n.begin(), verts_n.end(), verts_o.begin() );
00581         }
00582     }
00583 
00584     // write the tag back onto vertex coordinates
00585     rval = mb->set_coords( verts, &verts_n[0] );
00586     RC;
00587 
00588     dbgprint( "-- Final iterate error = " << global_max << ".\n" );
00589     return MB_SUCCESS;
00590 }
00591 
00592 /*
00593   Standard Laplacian Smooth Filter
00594 
00595   Additional references: http://www.doc.ic.ac.uk/~gr409/thesis-MSc.pdf
00596 */
00597 ErrorCode laplacianFilter( Core* mb, moab::ParallelComm* pcomm, moab::Range& verts, int dim, Tag fixed,
00598                            std::vector< double >& verts_o, std::vector< double >& verts_n, bool use_updated )
00599 {
00600     ErrorCode rval;
00601     std::vector< int > fix_tag( verts.size() );
00602     double* data;
00603 
00604     memcpy( &verts_n[0], &verts_o[0], sizeof( double ) * verts_o.size() );
00605 
00606     if( use_updated )
00607         data = &verts_n[0];
00608     else
00609         data = &verts_o[0];
00610 
00611     // filter verts down to owned ones and get fixed tag for them
00612     Range owned_verts;
00613     if( pcomm->size() > 1 )
00614     {
00615 #ifdef MOAB_HAVE_MPI
00616         rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );
00617         if( rval != MB_SUCCESS ) return rval;
00618 #endif
00619     }
00620     else
00621         owned_verts = verts;
00622 
00623     rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );
00624     RC;
00625 
00626     int vindex = 0;
00627     for( Range::const_iterator vit = owned_verts.begin(); vit != owned_verts.end(); ++vit, vindex++ )
00628     {
00629         // if !fixed
00630         if( fix_tag[vindex] ) continue;
00631 
00632         const int index = verts.index( *vit ) * 3;
00633 
00634         moab::Range adjverts, adjelems;
00635         // Find the neighboring vertices (1-ring neighborhood)
00636         rval = mb->get_adjacencies( &( *vit ), 1, dim, false, adjelems );
00637         RC;
00638         rval = mb->get_connectivity( adjelems, adjverts );
00639         RC;
00640         adjverts.erase( *vit );
00641 
00642         const int nadjs = adjverts.size();
00643         if( nadjs )
00644         {
00645             double delta[3] = { 0.0, 0.0, 0.0 };
00646 
00647             // Add the vertices and divide by the number of vertices
00648             for( int j = 0; j < nadjs; ++j )
00649             {
00650                 const int joffset = verts.index( adjverts[j] ) * 3;
00651                 delta[0] += data[joffset + 0];
00652                 delta[1] += data[joffset + 1];
00653                 delta[2] += data[joffset + 2];
00654             }
00655 
00656             verts_n[index + 0] = delta[0] / nadjs;
00657             verts_n[index + 1] = delta[1] / nadjs;
00658             verts_n[index + 2] = delta[2] / nadjs;
00659         }
00660     }
00661 
00662     return MB_SUCCESS;
00663 }
00664 
00665 /*
00666   HC (Humphrey’s Classes) Smooth Algorithm - Reduces Shrinkage of Laplacian Smoother
00667 
00668   Link:
00669   http://informatikbuero.com/downloads/Improved_Laplacian_Smoothing_of_Noisy_Surface_Meshes.pdf
00670 
00671   Where sv - original points
00672       pv - previous points,
00673       alpha [0..1] influences previous points pv, e.g. 0
00674       beta  [0..1] e.g. > 0.5
00675 */
00676 ErrorCode hcFilter( Core* mb, moab::ParallelComm* pcomm, moab::Range& verts, int dim, Tag fixed,
00677                     std::vector< double >& verts_o, std::vector< double >& verts_n, double alpha, double beta )
00678 {
00679     ErrorCode rval;
00680     std::vector< double > verts_hc( verts_o.size() );
00681     std::vector< int > fix_tag( verts.size() );
00682 
00683     // Perform Laplacian Smooth
00684     rval = laplacianFilter( mb, pcomm, verts, dim, fixed, verts_o, verts_n );
00685     RC;
00686 
00687     // Compute Differences
00688     for( unsigned index = 0; index < verts_o.size(); ++index )
00689     {
00690         verts_hc[index] = verts_n[index] - ( alpha * verts_n[index] + ( 1.0 - alpha ) * verts_o[index] );
00691     }
00692 
00693     // filter verts down to owned ones and get fixed tag for them
00694     Range owned_verts;
00695 #ifdef MOAB_HAVE_MPI
00696     if( pcomm->size() > 1 )
00697     {
00698         rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );
00699         if( rval != MB_SUCCESS ) return rval;
00700     }
00701     else
00702 #endif
00703         owned_verts = verts;
00704 
00705     rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );
00706     RC;
00707 
00708     int vindex = 0;
00709     for( Range::const_iterator vit = owned_verts.begin(); vit != owned_verts.end(); ++vit, vindex++ )
00710     {
00711         // if !fixed
00712         if( fix_tag[vindex] ) continue;
00713 
00714         const int index = verts.index( *vit ) * 3;
00715 
00716         moab::Range adjverts, adjelems;
00717         // Find the neighboring vertices (1-ring neighborhood)
00718         rval = mb->get_adjacencies( &( *vit ), 1, dim, false, adjelems );
00719         RC;
00720         rval = mb->get_connectivity( adjelems, adjverts );
00721         RC;
00722         adjverts.erase( *vit );
00723 
00724         const int nadjs = adjverts.size();
00725         double delta[3] = { 0.0, 0.0, 0.0 };
00726 
00727         for( int j = 0; j < nadjs; ++j )
00728         {
00729             const int joffset = verts.index( adjverts[j] ) * 3;
00730             delta[0] += verts_hc[joffset + 0];
00731             delta[1] += verts_hc[joffset + 1];
00732             delta[2] += verts_hc[joffset + 2];
00733         }
00734 
00735         verts_n[index + 0] -= beta * verts_hc[index + 0] + ( ( 1.0 - beta ) / nadjs ) * delta[0];
00736         verts_n[index + 1] -= beta * verts_hc[index + 1] + ( ( 1.0 - beta ) / nadjs ) * delta[1];
00737         verts_n[index + 2] -= beta * verts_hc[index + 2] + ( ( 1.0 - beta ) / nadjs ) * delta[2];
00738     }
00739 
00740     return MB_SUCCESS;
00741 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines