MOAB: Mesh Oriented datABase  (version 5.2.1)
LloydRelaxation.cpp
Go to the documentation of this file.
00001 /** @example LloydRelaxation.cpp \n
00002  * \brief Perform Lloyd relaxation on a mesh and its dual \n
00003  * <b>To run</b>: mpiexec -np <np> LloydRelaxation [filename]\n
00004  *
00005  * Briefly, Lloyd 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  * 'lloydfinal.h5m' in the current directory (H5M format must be used since the file is written in
00017  * parallel).
00018  */
00019 
00020 #include "moab/Core.hpp"
00021 #ifdef MOAB_HAVE_MPI
00022 #include "moab/ParallelComm.hpp"
00023 #include "MBParallelConventions.h"
00024 #endif
00025 #include "moab/Skinner.hpp"
00026 #include "moab/CN.hpp"
00027 #include "moab/CartVect.hpp"
00028 #include <iostream>
00029 #include <sstream>
00030 
00031 using namespace moab;
00032 using namespace std;
00033 
00034 string test_file_name = string( MESH_DIR ) + string( "/surfrandomtris-4part.h5m" );
00035 
00036 ErrorCode perform_lloyd_relaxation( Interface* mb, Range& verts, Range& cells, Tag fixed, int num_its, int report_its );
00037 
00038 int main( int argc, char** argv )
00039 {
00040     int num_its    = 10;
00041     int report_its = 1;
00042 
00043 #ifdef MOAB_HAVE_MPI
00044     MPI_Init( &argc, &argv );
00045 #endif
00046 
00047     // Need option handling here for input filename
00048     if( argc > 1 )
00049     {
00050         // User has input a mesh file
00051         test_file_name = argv[1];
00052     }
00053 
00054     // Get MOAB instance
00055     Interface* mb = new( std::nothrow ) Core;
00056     if( NULL == mb ) return 1;
00057 
00058 #ifdef MOAB_HAVE_MPI
00059     // Get the ParallelComm instance
00060     ParallelComm* pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
00061     int nprocs          = pcomm->size();
00062 #else
00063     int nprocs = 1;
00064 #endif
00065     string options;
00066     if( nprocs > 1 )  // If reading in parallel, need to tell it how
00067         options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
00068                   "PARALLEL_GHOSTS=2.0.1;DEBUG_IO=0;DEBUG_PIO=0";
00069 
00070     // Load the test file with specified options
00071     ErrorCode rval = mb->load_file( test_file_name.c_str(), 0, options.c_str() );MB_CHK_ERR( rval );
00072 
00073     // Make tag to specify fixed vertices, since it's input to the algorithm; use a default value of
00074     // non-fixed so we only need to set the fixed tag for skin vertices
00075     Tag fixed;
00076     int def_val = 0;
00077     rval        = mb->tag_get_handle( "fixed", 1, MB_TYPE_INTEGER, fixed, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );MB_CHK_ERR( rval );
00078 
00079     // Get all vertices and faces
00080     Range verts, faces, skin_verts;
00081     rval = mb->get_entities_by_type( 0, MBVERTEX, verts );MB_CHK_ERR( rval );
00082     rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );
00083 
00084     // Get the skin vertices of those faces and mark them as fixed; we don't want to fix the
00085     // vertices on a part boundary, but since we exchanged a layer of ghost faces, those vertices
00086     // aren't on the skin locally ok to mark non-owned skin vertices too, I won't move those anyway
00087     // use MOAB's skinner class to find the skin
00088     Skinner skinner( mb );
00089     rval = skinner.find_skin( 0, faces, true, skin_verts );MB_CHK_ERR( rval );  // 'true' param indicates we want vertices back, not faces
00090 
00091     vector< int > fix_tag( skin_verts.size(), 1 );  // Initialized to 1 to indicate fixed
00092     rval = mb->tag_set_data( fixed, skin_verts, &fix_tag[0] );MB_CHK_ERR( rval );
00093 
00094     // Now perform the Lloyd relaxation
00095     rval = perform_lloyd_relaxation( mb, verts, faces, fixed, num_its, report_its );MB_CHK_ERR( rval );
00096 
00097     // Delete fixed tag, since we created it here
00098     rval = mb->tag_delete( fixed );MB_CHK_ERR( rval );
00099 
00100     // Output file, using parallel write
00101 
00102 #ifdef MOAB_HAVE_MPI
00103     options = "PARALLEL=WRITE_PART";
00104 #endif
00105 
00106     rval = mb->write_file( "lloydfinal.h5m", NULL, options.c_str() );MB_CHK_ERR( rval );
00107 
00108     // Delete MOAB instance
00109     delete mb;
00110 
00111 #ifdef MOAB_HAVE_MPI
00112     MPI_Finalize();
00113 #endif
00114 
00115     return 0;
00116 }
00117 
00118 ErrorCode perform_lloyd_relaxation( Interface* mb, Range& verts, Range& faces, Tag fixed, int num_its, int report_its )
00119 {
00120     ErrorCode rval;
00121 
00122 #ifdef MOAB_HAVE_MPI
00123     ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
00124     int nprocs          = pcomm->size();
00125 #else
00126     int nprocs = 1;
00127 #endif
00128 
00129     // Perform Lloyd relaxation:
00130     // 1. Setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
00131 
00132     // Get all verts coords into tag; don't need to worry about filtering out fixed verts,
00133     // we'll just be setting to their fixed coords
00134     vector< double > vcentroids( 3 * verts.size() );
00135     rval = mb->get_coords( verts, &vcentroids[0] );MB_CHK_ERR( rval );
00136 
00137     Tag centroid;
00138     rval = mb->tag_get_handle( "centroid", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_ERR( rval );
00139     rval = mb->tag_set_data( centroid, verts, &vcentroids[0] );MB_CHK_ERR( rval );
00140 
00141     // Filter verts down to owned ones and get fixed tag for them
00142     Range owned_verts, shared_owned_verts;
00143     if( nprocs > 1 )
00144     {
00145 #ifdef MOAB_HAVE_MPI
00146         rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );MB_CHK_ERR( rval );
00147 #endif
00148     }
00149     else
00150         owned_verts = verts;
00151     vector< int > fix_tag( owned_verts.size() );
00152     rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );MB_CHK_ERR( rval );
00153 
00154     // Now fill vcentroids array with positions of just owned vertices, since those are the ones
00155     // we're actually computing
00156     vcentroids.resize( 3 * owned_verts.size() );
00157     rval = mb->tag_get_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );
00158 
00159 #ifdef MOAB_HAVE_MPI
00160     // Get shared owned verts, for exchanging tags
00161     rval = pcomm->get_shared_entities( -1, shared_owned_verts, 0, false, true );MB_CHK_ERR( rval );
00162     // Workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging
00163     // tags for all shared entities
00164     if( shared_owned_verts.empty() ) shared_owned_verts.insert( *verts.begin() );
00165 #endif
00166 
00167     // Some declarations for later iterations
00168     vector< double > fcentroids( 3 * faces.size() );         // fcentroids for face centroids
00169     vector< double > ctag( 3 * CN::MAX_NODES_PER_ELEMENT );  // Temporary coordinate storage for verts bounding a face
00170     const EntityHandle* conn;                                // const ptr & size to face connectivity
00171     int nconn;
00172     Range::iterator fit, vit;          // For iterating over faces, verts
00173     int f, v;                          // For indexing into centroid vectors
00174     vector< EntityHandle > adj_faces;  // Used in vertex iteration
00175 
00176     // 2. For num_its iterations:
00177     for( int nit = 0; nit < num_its; nit++ )
00178     {
00179         double mxdelta = 0.0;
00180 
00181         // 2a. For each face: centroid = sum(vertex centroids)/num_verts_in_cell
00182         for( fit = faces.begin(), f = 0; fit != faces.end(); ++fit, f++ )
00183         {
00184             // Get verts for this face
00185             rval = mb->get_connectivity( *fit, conn, nconn );MB_CHK_ERR( rval );
00186             // Get centroid tags for those verts
00187             rval = mb->tag_get_data( centroid, conn, nconn, &ctag[0] );MB_CHK_ERR( rval );
00188             fcentroids[3 * f + 0] = fcentroids[3 * f + 1] = fcentroids[3 * f + 2] = 0.0;
00189             for( v = 0; v < nconn; v++ )
00190             {
00191                 fcentroids[3 * f + 0] += ctag[3 * v + 0];
00192                 fcentroids[3 * f + 1] += ctag[3 * v + 1];
00193                 fcentroids[3 * f + 2] += ctag[3 * v + 2];
00194             }
00195             for( v = 0; v < 3; v++ )
00196                 fcentroids[3 * f + v] /= nconn;
00197         }
00198         rval = mb->tag_set_data( centroid, faces, &fcentroids[0] );MB_CHK_ERR( rval );
00199 
00200         // 2b. For each owned vertex:
00201         for( vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); ++vit, v++ )
00202         {
00203             // if !fixed
00204             if( fix_tag[v] ) continue;
00205             // vertex centroid = sum(cell centroids)/ncells
00206             adj_faces.clear();
00207             rval = mb->get_adjacencies( &( *vit ), 1, 2, false, adj_faces );MB_CHK_ERR( rval );
00208             rval = mb->tag_get_data( centroid, &adj_faces[0], adj_faces.size(), &fcentroids[0] );MB_CHK_ERR( rval );
00209             double vnew[] = { 0.0, 0.0, 0.0 };
00210             for( f = 0; f < (int)adj_faces.size(); f++ )
00211             {
00212                 vnew[0] += fcentroids[3 * f + 0];
00213                 vnew[1] += fcentroids[3 * f + 1];
00214                 vnew[2] += fcentroids[3 * f + 2];
00215             }
00216             for( f = 0; f < 3; f++ )
00217                 vnew[f] /= adj_faces.size();
00218             double delta = ( CartVect( vnew ) - CartVect( &vcentroids[3 * v] ) ).length();
00219             mxdelta      = std::max( delta, mxdelta );
00220             for( f = 0; f < 3; f++ )
00221                 vcentroids[3 * v + f] = vnew[f];
00222         }
00223 
00224         // Set the centroid tag; having them only in vcentroids array isn't enough, as vertex
00225         // centroids are accessed randomly in loop over faces
00226         rval = mb->tag_set_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );
00227 
00228         // 2c. Exchange tags on owned verts
00229         if( nprocs > 1 )
00230         {
00231 #ifdef MOAB_HAVE_MPI
00232             rval = pcomm->exchange_tags( centroid, shared_owned_verts );MB_CHK_ERR( rval );
00233 #endif
00234         }
00235 
00236         if( !( nit % report_its ) )
00237         {
00238             // Global reduce for maximum delta, then report it
00239             double global_max = mxdelta;
00240             int myrank        = 0;
00241 #ifdef MOAB_HAVE_MPI
00242             if( nprocs > 1 ) MPI_Reduce( &mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
00243             myrank = pcomm->rank();
00244 #endif
00245             if( 1 == nprocs || 0 == myrank ) cout << "Max delta = " << global_max << endl;
00246         }
00247     }
00248 
00249     // Write the tag back onto vertex coordinates
00250     rval = mb->set_coords( owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );
00251 
00252     // Delete the centroid tag, since we don't need it anymore
00253     rval = mb->tag_delete( centroid );MB_CHK_ERR( rval );
00254 
00255     return MB_SUCCESS;
00256 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines