LloydRelaxation.cpp File Reference

Example demonstrating Lloyd relaxation for mesh smoothing. More...

#include "moab/Core.hpp"
#include "moab/Skinner.hpp"
#include "moab/CN.hpp"
#include "moab/CartVect.hpp"
#include <iostream>
#include <sstream>

Include dependency graph for LloydRelaxation.cpp:

Go to the source code of this file.

Functions
ErrorCode	perform_lloyd_relaxation (Interface *mb, Range &verts, Range &cells, Tag fixed, int num_its, int report_its)
int	main (int argc, char **argv)

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

Detailed Description

Example demonstrating Lloyd relaxation for mesh smoothing.

This example shows how to:

• Load a parallel mesh with ghost layers
• Perform Lloyd relaxation iterations for mesh smoothing
• Handle fixed vertices (e.g., boundary vertices)
• Exchange ghost data between processors
• Measure convergence based on vertex movement
• Write the smoothed mesh to a parallel file

Lloyd relaxation is a technique to smooth out a mesh by iteratively moving vertices to the centroids of their connected cells.

Author

MOAB Development Team

Date

2024

Perform Lloyd relaxation on a mesh and its dual
To run: mpiexec -np np LloydRelaxation [filename]
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).

Parameters

argc	Number of command line arguments
argv	Command line arguments array

Returns

0 on success, 1 on failure

Definition in file LloydRelaxation.cpp.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 60 of file LloydRelaxation.cpp.

{
    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
    MB_CHK_ERR( mb->load_file( test_file_name.c_str(), 0, options.c_str() ) );
 
    // 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;
    MB_CHK_ERR( mb->tag_get_handle( "fixed", 1, MB_TYPE_INTEGER, fixed, MB_TAG_CREAT | MB_TAG_DENSE, &def_val ) );
 
    // Get all vertices and faces
    Range verts, faces, skin_verts;
    MB_CHK_ERR( mb->get_entities_by_type( 0, MBVERTEX, verts ) );
    MB_CHK_ERR( mb->get_entities_by_dimension( 0, 2, faces ) );
 
    // 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 );
    MB_CHK_ERR(
        skinner.find_skin( 0, faces, true, skin_verts ) );  // 'true' param indicates we want vertices back, not faces
 
    vector< int > fix_tag( skin_verts.size(), 1 );  // Initialized to 1 to indicate fixed
    MB_CHK_ERR( mb->tag_set_data( fixed, skin_verts, &fix_tag[0] ) );
 
    // Now perform the Lloyd relaxation
    MB_CHK_ERR( perform_lloyd_relaxation( mb, verts, faces, fixed, num_its, report_its ) );
 
    // Delete fixed tag, since we created it here
    MB_CHK_ERR( mb->tag_delete( fixed ) );
 
    // Output file, using parallel write
 
#ifdef MOAB_HAVE_MPI
    options = "PARALLEL=WRITE_PART";
#endif
 
    MB_CHK_ERR( mb->write_file( "lloydfinal.h5m", NULL, options.c_str() ) );
 
    // Delete MOAB instance
    delete mb;
 
#ifdef MOAB_HAVE_MPI
    MPI_Finalize();
#endif
 
    return 0;
}

References moab::Skinner::find_skin(), moab::Core::get_entities_by_dimension(), moab::Core::get_entities_by_type(), moab::Core::load_file(), mb, MB_CHK_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_INTEGER, MBVERTEX, perform_lloyd_relaxation(), moab::Range::size(), moab::ParallelComm::size(), moab::Core::tag_delete(), moab::Core::tag_get_handle(), moab::Core::tag_set_data(), test_file_name, and moab::Core::write_file().

perform_lloyd_relaxation()

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

Definition at line 141 of file LloydRelaxation.cpp.

{
    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() );
    MB_CHK_ERR( mb->get_coords( verts, &vcentroids[0] ) );
 
    Tag centroid;
    MB_CHK_ERR( mb->tag_get_handle( "centroid", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE ) );
    MB_CHK_ERR( mb->tag_set_data( centroid, verts, &vcentroids[0] ) );
 
    // 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
        MB_CHK_ERR( pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts ) );
#endif
    }
    else
        owned_verts = verts;
    vector< int > fix_tag( owned_verts.size() );
    MB_CHK_ERR( mb->tag_get_data( fixed, owned_verts, &fix_tag[0] ) );
 
    // 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() );
    MB_CHK_ERR( mb->tag_get_data( centroid, owned_verts, &vcentroids[0] ) );
 
#ifdef MOAB_HAVE_MPI
    // Get shared owned verts, for exchanging tags
    MB_CHK_ERR( pcomm->get_shared_entities( -1, shared_owned_verts, 0, false, true ) );
    // 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
            MB_CHK_ERR( mb->get_connectivity( *fit, conn, nconn ) );
            // Get centroid tags for those verts
            MB_CHK_ERR( mb->tag_get_data( centroid, conn, nconn, &ctag[0] ) );
            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;
        }
        MB_CHK_ERR( mb->tag_set_data( centroid, faces, &fcentroids[0] ) );
 
        // 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();
            MB_CHK_ERR( mb->get_adjacencies( &( *vit ), 1, 2, false, adj_faces ) );
            MB_CHK_ERR( mb->tag_get_data( centroid, &adj_faces[0], adj_faces.size(), &fcentroids[0] ) );
            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
        MB_CHK_ERR( mb->tag_set_data( centroid, owned_verts, &vcentroids[0] ) );
 
        // 2c. Exchange tags on owned verts
        if( nprocs > 1 )
        {
#ifdef MOAB_HAVE_MPI
            MB_CHK_ERR( pcomm->exchange_tags( centroid, shared_owned_verts ) );
#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
    MB_CHK_ERR( mb->set_coords( owned_verts, &vcentroids[0] ) );
 
    // Delete the centroid tag, since we don't need it anymore
    MB_CHK_ERR( mb->tag_delete( centroid ) );
 
    return MB_SUCCESS;
}

References moab::Range::begin(), moab::ParallelComm::comm(), moab::Range::empty(), moab::Range::end(), ErrorCode, moab::ParallelComm::exchange_tags(), moab::ParallelComm::filter_pstatus(), moab::Core::get_adjacencies(), moab::Core::get_connectivity(), moab::Core::get_coords(), moab::ParallelComm::get_pcomm(), moab::ParallelComm::get_shared_entities(), moab::Range::insert(), length(), moab::CN::MAX_NODES_PER_ELEMENT, mb, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, PSTATUS_NOT, PSTATUS_NOT_OWNED, moab::ParallelComm::rank(), moab::Core::set_coords(), moab::Range::size(), moab::ParallelComm::size(), moab::Core::tag_delete(), moab::Core::tag_get_data(), moab::Core::tag_get_handle(), and moab::Core::tag_set_data().

Referenced by main().

Variable Documentation

test_file_name

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

Definition at line 56 of file LloydRelaxation.cpp.

Referenced by main().