LloydRelaxation.cpp

Go to the documentation of this file.

/** @example LloydRelaxation.cpp \n
 * \brief Perform Lloyd relaxation on a mesh and its dual \n
 * <b>To run</b>: mpiexec -np <np> LloydRelaxation [filename]\n
 *
 * Briefly, Lloyd 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
 * 'lloydfinal.h5m' in the current directory (H5M format must be used since the file is written in
 * parallel).
 */

#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/CartVect.hpp"
#include <iostream>
#include <sstream>

using namespace moab;
using namespace std;

string test_file_name = string( MESH_DIR ) + string( "/surfrandomtris-4part.h5m" );

ErrorCode perform_lloyd_relaxation( Interface* mb, Range& verts, Range& cells, Tag fixed, int num_its, int report_its );

int main( int argc, char** argv )
{
    int num_its    = 10;
    int report_its = 1;

#ifdef MOAB_HAVE_MPI
    MPI_Init( &argc, &argv );
#endif

    // Need option handling here for input filename
    if( argc > 1 )
    {
        // User has input a mesh file
        test_file_name = argv[1];
    }

    // Get MOAB instance
    Interface* mb = new( std::nothrow ) Core;
    if( NULL == mb ) return 1;

#ifdef MOAB_HAVE_MPI
    // Get the ParallelComm instance
    ParallelComm* pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
    int nprocs          = pcomm->size();
#else
    int nprocs = 1;
#endif
    string options;
    if( nprocs > 1 )  // If reading in parallel, need to tell it how
        options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
                  "PARALLEL_GHOSTS=2.0.1;DEBUG_IO=0;DEBUG_PIO=0";

    // Load the test file with specified options
    ErrorCode rval = mb->load_file( test_file_name.c_str(), 0, options.c_str() );MB_CHK_ERR( rval );

    // 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 );MB_CHK_ERR( rval );

    // Get all vertices and faces
    Range verts, faces, skin_verts;
    rval = mb->get_entities_by_type( 0, MBVERTEX, verts );MB_CHK_ERR( rval );
    rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );

    // Get the skin vertices of those faces 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 faces, 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( 0, faces, true, skin_verts );MB_CHK_ERR( rval );  // 'true' param indicates we want vertices back, not faces

    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] );MB_CHK_ERR( rval );

    // Now perform the Lloyd relaxation
    rval = perform_lloyd_relaxation( mb, verts, faces, fixed, num_its, report_its );MB_CHK_ERR( rval );

    // Delete fixed tag, since we created it here
    rval = mb->tag_delete( fixed );MB_CHK_ERR( rval );

    // Output file, using parallel write

#ifdef MOAB_HAVE_MPI
    options = "PARALLEL=WRITE_PART";
#endif

    rval = mb->write_file( "lloydfinal.h5m", NULL, options.c_str() );MB_CHK_ERR( rval );

    // Delete MOAB instance
    delete mb;

#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif

    return 0;
}

ErrorCode perform_lloyd_relaxation( Interface* mb, Range& verts, Range& faces, Tag fixed, int num_its, int report_its )
{
    ErrorCode rval;

#ifdef MOAB_HAVE_MPI
    ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
    int nprocs          = pcomm->size();
#else
    int nprocs = 1;
#endif

    // Perform Lloyd 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
    vector< double > vcentroids( 3 * verts.size() );
    rval = mb->get_coords( verts, &vcentroids[0] );MB_CHK_ERR( rval );

    Tag centroid;
    rval = mb->tag_get_handle( "centroid", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_ERR( rval );
    rval = mb->tag_set_data( centroid, verts, &vcentroids[0] );MB_CHK_ERR( rval );

    // Filter verts down to owned ones and get fixed tag for them
    Range owned_verts, shared_owned_verts;
    if( nprocs > 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;
    vector< int > fix_tag( owned_verts.size() );
    rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );MB_CHK_ERR( rval );

    // Now fill vcentroids array with positions of just owned vertices, since those are the ones
    // we're actually computing
    vcentroids.resize( 3 * owned_verts.size() );
    rval = mb->tag_get_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );

#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

    // Some declarations for later iterations
    vector< double > fcentroids( 3 * faces.size() );         // fcentroids for face centroids
    vector< double > ctag( 3 * CN::MAX_NODES_PER_ELEMENT );  // Temporary coordinate storage for verts bounding a face
    const EntityHandle* conn;                                // const ptr & size to face connectivity
    int nconn;
    Range::iterator fit, vit;          // For iterating over faces, verts
    int f, v;                          // For indexing into centroid vectors
    vector< EntityHandle > adj_faces;  // Used in vertex iteration

    // 2. For num_its iterations:
    for( int nit = 0; nit < num_its; nit++ )
    {
        double mxdelta = 0.0;

        // 2a. For each face: centroid = sum(vertex centroids)/num_verts_in_cell
        for( fit = faces.begin(), f = 0; fit != faces.end(); ++fit, f++ )
        {
            // Get verts for this face
            rval = mb->get_connectivity( *fit, conn, nconn );MB_CHK_ERR( rval );
            // Get centroid tags for those verts
            rval = mb->tag_get_data( centroid, conn, nconn, &ctag[0] );MB_CHK_ERR( rval );
            fcentroids[3 * f + 0] = fcentroids[3 * f + 1] = fcentroids[3 * f + 2] = 0.0;
            for( v = 0; v < nconn; v++ )
            {
                fcentroids[3 * f + 0] += ctag[3 * v + 0];
                fcentroids[3 * f + 1] += ctag[3 * v + 1];
                fcentroids[3 * f + 2] += ctag[3 * v + 2];
            }
            for( v = 0; v < 3; v++ )
                fcentroids[3 * f + v] /= nconn;
        }
        rval = mb->tag_set_data( centroid, faces, &fcentroids[0] );MB_CHK_ERR( rval );

        // 2b. For each owned vertex:
        for( vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); ++vit, v++ )
        {
            // if !fixed
            if( fix_tag[v] ) continue;
            // vertex centroid = sum(cell centroids)/ncells
            adj_faces.clear();
            rval = mb->get_adjacencies( &( *vit ), 1, 2, false, adj_faces );MB_CHK_ERR( rval );
            rval = mb->tag_get_data( centroid, &adj_faces[0], adj_faces.size(), &fcentroids[0] );MB_CHK_ERR( rval );
            double vnew[] = { 0.0, 0.0, 0.0 };
            for( f = 0; f < (int)adj_faces.size(); f++ )
            {
                vnew[0] += fcentroids[3 * f + 0];
                vnew[1] += fcentroids[3 * f + 1];
                vnew[2] += fcentroids[3 * f + 2];
            }
            for( f = 0; f < 3; f++ )
                vnew[f] /= adj_faces.size();
            double delta = ( CartVect( vnew ) - CartVect( &vcentroids[3 * v] ) ).length();
            mxdelta      = std::max( delta, mxdelta );
            for( f = 0; f < 3; f++ )
                vcentroids[3 * v + f] = vnew[f];
        }

        // Set the centroid tag; having them only in vcentroids array isn't enough, as vertex
        // centroids are accessed randomly in loop over faces
        rval = mb->tag_set_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );

        // 2c. Exchange tags on owned verts
        if( nprocs > 1 )
        {
#ifdef MOAB_HAVE_MPI
            rval = pcomm->exchange_tags( centroid, shared_owned_verts );MB_CHK_ERR( rval );
#endif
        }

        if( !( nit % report_its ) )
        {
            // Global reduce for maximum delta, then report it
            double global_max = mxdelta;
            int myrank        = 0;
#ifdef MOAB_HAVE_MPI
            if( nprocs > 1 ) MPI_Reduce( &mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
            myrank = pcomm->rank();
#endif
            if( 1 == nprocs || 0 == myrank ) cout << "Max delta = " << global_max << endl;
        }
    }

    // Write the tag back onto vertex coordinates
    rval = mb->set_coords( owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );

    // Delete the centroid tag, since we don't need it anymore
    rval = mb->tag_delete( centroid );MB_CHK_ERR( rval );

    return MB_SUCCESS;
}