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