ExchangeHalos.cpp

Example program that shows the use case for performing tag data exchange between parallel processors in order to sync data on shared entities.

This example :

1. Initialize MPI and instantiates MOAB
2. Gets user options: Input mesh file name, vector tag length, ghost layer size etc
3. Create the root and partition sets
4. Instantiate ParallelComm and read the mesh file in parallel using appropriate options
5. Create the required number of ghost layers as requested by the user (default = 3)
6. Get 2D MPAS polygonal entities in the mesh and filter to get only the "owned" entities
7. Create two tags: scalar_variable (single data/cell) and vector_variable (multiple data/cell)
8. Set tag data using analytical functions for both scalar and vector fields on owned entities
9. Exchange shared entity information and tags between processors
   1. If debugging is turned on, store mesh file and tag on root process (will not contain data on shared entities)
   2. Perform exchange of scalar tag data on shared entities
   3. Perform exchange of vector tag data on shared entities
   4. If debugging is turned on, store mesh file and tag on root process (will now contain data on all entities)
10. Destroy the MOAB instance and finalize MPI

To run: mpiexec -n np ./ExchangeHalos –input mpas_mesh_file –nghosts ghostlayers –vtaglength vector_component_size –nexchanges number_of_exchange_runs Example: mpiexec -n 4 ./ExchangeHalos –input $MOAB_DIR/MeshFiles/unittest/io/mpasx1.642.t.2.nc –nghosts 3 –vtaglength 100

NOTE: –debug option can be added to write out extra files in h5m format to visualize some output (written from root task only)

/** @example ExchangeHalos.cpp
 * 
 * Example program that shows the use case for performing tag data exchange
 * between parallel processors in order to sync data on shared entities. \n\n
 *
 * <b>This example </b>:
 *    -# Initialize MPI and instantiates MOAB
 *    -# Gets user options: Input mesh file name, vector tag length, ghost layer size etc
 *    -# Create the root and partition sets
 *    -# Instantiate ParallelComm and read the mesh file in parallel using appropriate options
 *    -# Create the required number of ghost layers as requested by the user (default = 3)
 *    -# Get 2D MPAS polygonal entities in the mesh and filter to get only the "owned" entities
 *    -# Create two tags: scalar_variable (single data/cell) and vector_variable (multiple data/cell)
 *    -# Set tag data using analytical functions for both scalar and vector fields on owned entities
 *    -# Exchange shared entity information and tags between processors
 *      -# If debugging is turned on, store mesh file and tag on root process (will not contain data on shared entities)
 *      -# Perform exchange of scalar tag data on shared entities
 *      -# Perform exchange of vector tag data on shared entities
 *      -# If debugging is turned on, store mesh file and tag on root process (will now contain data on *all* entities)
 *    -#  Destroy the MOAB instance and finalize MPI
 *
 * <b>To run: </b>
 *      mpiexec -n np ./ExchangeHalos --input \c mpas_mesh_file --nghosts \c ghostlayers --vtaglength \c vector_component_size 
 *                    --nexchanges \c number_of_exchange_runs
 * <b>Example:</b>
 *      mpiexec -n 4 ./ExchangeHalos --input $MOAB_DIR/MeshFiles/unittest/io/mpasx1.642.t.2.nc --nghosts 3 --vtaglength 100
 *
 * NOTE: --debug option can be added to write out extra files in h5m format to visualize some output (written from root
 * task only)
 *
 */
 
// MOAB includes
#include "moab/Core.hpp"
#include "moab/CpuTimer.hpp"
#include "moab/ProgOptions.hpp"
 
#ifndef MOAB_HAVE_MPI
#error "Please build MOAB with MPI..."
#endif
 
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
 
// C++ includes
#include <iostream>
#include <string>
 
using namespace moab;
using namespace std;
#define dbgprint( MSG )                                           \
    do                                                            \
    {                                                             \
        if( context.proc_id == 0 ) std::cout << MSG << std::endl; \
    } while( false )
 
#define runchk( CODE, MSG )         \
    do                              \
    {                               \
        moab::ErrorCode err = CODE; \
        MB_CHK_SET_ERR( err, MSG ); \
    } while( false )
 
#define runchk_cont( CODE, MSG )                               \
    do                                                         \
    {                                                          \
        moab::ErrorCode err = CODE;                            \
        MB_CHK_ERR_CONT( err );                                \
        if( err ) std::cout << "Error:: " << MSG << std::endl; \
    } while( false )
 
// Run time context structure
 
/// @brief The RunttimeContext is an example specific class to store
/// the run specific input data, MOAB datastructures used during the run
/// and provides other utility functions to profile operations etc
struct RuntimeContext
{
  public:
    int dimension{ 2 };           /// dimension of the problem
    std::string input_filename;   /// input file name (nc format)
    std::string output_filename;  /// output file name (h5m format)
    int ghost_layers{ 3 };        /// number of ghost layers
    std::string scalar_tagname;   /// scalar tag name
    std::string vector_tagname;   /// vector tag name
    int vector_length{ 3 };       /// length of the vector tag components
    int num_max_exchange{ 10 };   /// total number of exchange iterations
    bool debug_output{ false };   /// write debug output information?
    int proc_id{ 1 };             /// process identifier
    int num_procs{ 1 };           /// total number of processes
    double last_counter{ 0.0 };   /// last time counter between push/pop timer
 
    // MOAB objects
    moab::Core moab_interface;
    moab::ParallelComm* parallel_communicator;
    moab::EntityHandle fileset{ 0 }, partnset{ 0 };
 
    /// @brief Constructor: allocate MOAB interface and communicator, and initialize
    /// other data members with some default values
    explicit RuntimeContext( MPI_Comm comm = MPI_COMM_WORLD );
 
    /// @brief Destructor: deallocate MOAB interface and communicator
    ~RuntimeContext();
 
    /// @brief Parse the runtime command line options
    /// @param argc - number of command line arguments
    /// @param argv - command line arguments as string list
    void ParseCLOptions( int argc, char* argv[] );
 
    /// @brief Measure and start the timer to profile a task
    /// @param operation String name of the task being measured
    inline void timer_push( const std::string& operation );
 
    /// @brief Stop the timer and store the elapsed duration
    /// @param nruns Optional argument used to average the measured time
    void timer_pop( const int nruns = 1 );
 
    /// @brief Return the last elapsed time
    /// @return last_counter from timer_pop was called
    inline double last_elapsed() const;
 
    /// @brief Load a MOAB supported file (h5m or nc format) from disk
    ///        representing an MPAS mesh
    /// @param load_ghosts Optional boolean to specify whether to load ghosts
    ///                    when reading the file (only relevant for h5m)
    /// @return Error code if any (else MB_SUCCESS)
    moab::ErrorCode load_file( bool load_ghosts = false );
 
    /// @brief Create scalar and vector tags in the MOAB mesh instance
    /// @param tagScalar Tag reference to the scalar field
    /// @param tagVector Tag reference to the vector field
    /// @param entities Entities on which both the scalar and vector fields are defined
    /// @return Error code if any (else MB_SUCCESS)
    moab::ErrorCode create_sv_tags( moab::Tag& tagScalar, moab::Tag& tagVector, moab::Range& entities );
 
    /// @brief Evaluate some closed-form Spherical Harmonic functions with an optional multiplier term
    /// @param lon Longitude in lat-lon space
    /// @param lat Latitude  in lat-lon space
    /// @param type Function type
    /// @param multiplier Optional multiplier to scale value (default=1.0)
    /// @return value of the evaluated function
    inline double evaluate_function( double lon, double lat, int type = 1, double multiplier = 1.0 ) const
    {
        switch( type )
        {
            case 1:
                return ( 2.0 + std::pow( sin( 2.0 * lat ), 16.0 ) * cos( 16.0 * lon ) ) * multiplier;
            default:
                return ( 2.0 + cos( lon ) * cos( lon ) * cos( 2.0 * lat ) ) * multiplier;
        }
    }
 
  private:
    /// @brief Compute the centroids of elements in 2D lat/lon space
    /// @param entities Entities to compute centroids
    /// @return Vector of centroids (as lat/lon)
    std::vector< double > compute_centroids( const moab::Range& entities ) const;
 
    moab::CpuTimer mTimer;
    double mTimerOps{ 0.0 };
    std::string mOpName;
};
 
//
// Start of main test program
//
int main( int argc, char** argv )
{
    // Initialize MPI first
    MPI_Init( &argc, &argv );
 
    {
        // Create our context for this example run
        RuntimeContext context;
        dbgprint( "********** Exchange halos example **********\n" );
 
        // Get the input options
        context.ParseCLOptions( argc, argv );
 
        /////////////////////////////////////////////////////////////////////////
        // Print out the input parameters in use
        dbgprint( " -- Input Parameters -- " );
        dbgprint( "    Number of Processes  = " << context.num_procs );
        dbgprint( "    Input mesh           = " << context.input_filename );
        dbgprint( "    Ghost Layers         = " << context.ghost_layers );
        dbgprint( "    Scalar Tag name      = " << context.scalar_tagname );
        dbgprint( "    Vector Tag name      = " << context.vector_tagname );
        dbgprint( "    Vector Tag length    = " << context.vector_length << endl );
        /////////////////////////////////////////////////////////////////////////
 
        // Timer storage for all phases
        double elapsed_times[4];
 
        // Read the input file specified by user, in parallel, using appropriate options
        // Supports reading partitioned h5m files and MPAS nc files directly with online Zoltan partitioning
        context.timer_push( "Read input file" );
        {
            // Load the file from disk with given options
            runchk( context.load_file( false ), "MOAB::load_file failed for filename: " << context.input_filename );
        }
        context.timer_pop();
        elapsed_times[0] = context.last_elapsed();
 
        // Let the actual measurements begin...
        dbgprint( "\n- Starting execution -\n" );
 
        // We need to set up the ghost layers requested by the user. First correct for thin layers and then
        // call `exchange_ghost_cells` to prepare the mesh for use with halo regions
        context.timer_push( "Setup ghost layers" );
        {
            // Loop over the number of ghost layers needed and ask MOAB for layers 1 at a time
            for( int ighost = 0; ighost < context.ghost_layers; ++ighost )
            {
                // Exchange ghost cells
                int ghost_dimension  = context.dimension;
                int bridge_dimension = context.dimension - 1;
                // Let us now get all ghost layers from adjacent parts
                runchk( context.parallel_communicator->exchange_ghost_cells(
                            ghost_dimension, bridge_dimension, ( ighost + 1 ), 0, true /* store_remote_handles */,
                            true /* wait_all */, &context.fileset ),
                        "Exchange ghost cells failed" );  // true to store remote handles
 
                // Ensure that all processes understand about multi-shared vertices and entities
                // in case some adjacent parts are only m layers thick (where m < context.ghost_layers)
                if( ighost < context.ghost_layers - 1 )
                    runchk( context.parallel_communicator->correct_thin_ghost_layers(),
                            "Thin layer correction failed" );
            }
        }
        context.timer_pop();
        elapsed_times[1] = context.last_elapsed();
 
        // Get the 2D MPAS elements and filter it so that we have only owned elements
        Range dimEnts;
        {
            // Get all entities of dimension = dim
            runchk( context.moab_interface.get_entities_by_dimension( context.fileset, context.dimension, dimEnts ),
                    "Getting 2D entities failed" );
            // Get only owned entities! The ghosted/shared entities will get their data when we exchange
            // So let us filter entities based on the status: NOT x NOT_OWNED = OWNED status :-)
            runchk( context.parallel_communicator->filter_pstatus( dimEnts, PSTATUS_NOT_OWNED, PSTATUS_NOT ),
                    "Filtering pstatus failed" );
 
            // Aggregate the total number of elements in the mesh
            auto numEntities     = dimEnts.size();
            int numTotalEntities = 0;
            MPI_Reduce( &numEntities, &numTotalEntities, 1, MPI_INT, MPI_SUM, 0,
                        context.parallel_communicator->proc_config().proc_comm() );
 
            // We expect the total number of elements to be constant, immaterial of number of processes.
            // If not, we have a bug!
            dbgprint( "Total number of " << context.dimension << "D elements in the mesh = " << numTotalEntities );
        }
 
        Tag tagScalar = nullptr;
        Tag tagVector = nullptr;
        // Create two tag handles: scalar_variable and vector_variable
        // Set these tags with appropriate closed form functional data
        // based on element centroid information
        runchk( context.create_sv_tags( tagScalar, tagVector, dimEnts ), "Unable to create scalar and vector tags" );
 
        // let us write out the local mesh before tag_exchange is called
        // we expect to see data only on the owned entities - and ghosted entities should have default values
        if( context.debug_output && ( context.proc_id == 0 ) )  // only on root process, for debugging
        {
            dbgprint( "> Writing to file *before* ghost exchange " );
            runchk( context.moab_interface.write_file( "exchangeHalos_output_rank0_pre.h5m", "H5M", "" ),
                    "Writing to disk failed" );
        }
 
        // Perform exchange of tag data between neighboring tasks
        dbgprint( "> Exchanging tags between processors " );
        context.timer_push( "Exchange scalar tag data" );
        for( auto irun = 0; irun < context.num_max_exchange; ++irun )
        {
            // Exchange scalar tags between processors
            runchk( context.parallel_communicator->exchange_tags( tagScalar, dimEnts ),
                    "Exchanging scalar tag between processors failed" );
        }
        context.timer_pop( context.num_max_exchange );
        elapsed_times[2] = context.last_elapsed();
 
        context.timer_push( "Exchange vector tag data" );
        for( auto irun = 0; irun < context.num_max_exchange; ++irun )
        {
            // Exchange vector tags between processors
            runchk( context.parallel_communicator->exchange_tags( tagVector, dimEnts ),
                    "Exchanging vector tag between processors failed" );
        }
        context.timer_pop( context.num_max_exchange );
        elapsed_times[3] = context.last_elapsed();
 
        // let us write out the local mesh after tag_exchange is called
        // we expect to see real data on both owned and ghost entities in halo regions (non-default values)
        if( context.debug_output && ( context.proc_id == 0 ) )  // only on root process, for debugging
        {
            dbgprint( "> Writing to file *after* ghost exchange " );
            runchk( context.moab_interface.write_file( "exchangeHalos_output_rank0_post.h5m", "H5M", "" ),
                    "Writing to disk failed" );
        }
 
        // Write out the final mesh with the tag data and mesh -- just for verification
        if( context.debug_output )
        {
            dbgprint( "> Writing out the final mesh and data in MOAB h5m format. File = " << context.output_filename );
            string write_options = ( context.num_procs > 1 ? "PARALLEL=WRITE_PART;DEBUG_IO=0;" : "" );
            // Write out to output file to visualize reduction/exchange of tag data
            runchk( context.moab_interface.write_file( context.output_filename.c_str(), "H5M", write_options.c_str() ),
                    "File write failed" );
        }
 
        // Consolidated timing results: the data is listed as follows
        // [ntasks,  nghosts,  load_mesh(I/O),  exchange_ghost_cells(setup), exchange_tags(scalar),
        // exchange_tags(vector)]
        dbgprint( "\n> Consolidated: [" << context.num_procs << ", " << context.ghost_layers << ", " << elapsed_times[0]
                                        << ", " << elapsed_times[1] << ", " << elapsed_times[2] << ", "
                                        << elapsed_times[3] << "]," );
 
        // execution finished
        dbgprint( "\n********** ExchangeHalos Example DONE! **********" );
    }
    // Done, cleanup
    MPI_Finalize();
 
    return 0;
}
 
/// Implementation details ///
 
///
RuntimeContext::RuntimeContext( MPI_Comm comm )
    : input_filename( std::string( MESH_DIR ) + std::string( "/io/mpasx1.642.t.2.nc" ) ),
      output_filename( "exchangeHalos_output.h5m" ), scalar_tagname( "scalar_variable" ),
      vector_tagname( "vector_variable" )
{
    // Create sets for the mesh and partition.  Then pass these to the load_file functions to populate the mesh.
    runchk_cont( moab_interface.create_meshset( moab::MESHSET_SET, fileset ), "Creating root set failed" );
    runchk_cont( moab_interface.create_meshset( moab::MESHSET_SET, partnset ), "Creating partition set failed" );
 
    // Create the parallel communicator object with the partition handle associated with MOAB
    parallel_communicator = moab::ParallelComm::get_pcomm( &moab_interface, partnset, &comm );
 
    proc_id   = parallel_communicator->rank();
    num_procs = parallel_communicator->size();
}
 
RuntimeContext::~RuntimeContext()
{
    delete parallel_communicator;
}
 
void RuntimeContext::ParseCLOptions( int argc, char* argv[] )
{
    ProgOptions opts;
    // Input mesh
    opts.addOpt< std::string >( "input", "Input mesh filename to load in parallel. Default=data/default_mesh_holes.h5m",
                                &input_filename );
    // Output mesh
    opts.addOpt< void >( "debug", "Should we write output file? Default=false", &debug_output );
    opts.addOpt< std::string >( "output",
                                "Output mesh filename for verification (use --debug). Default=exchangeHalos_output.h5m",
                                &output_filename );
    // Dimension of the input mesh
    // Vector tag length
    opts.addOpt< int >( "vtaglength", "Size of vector components per each entity. Default=3", &vector_length );
    // Number of halo (ghost) regions
    opts.addOpt< int >( "nghosts", "Number of ghost layers (halos) to exchange. Default=3", &ghost_layers );
    // Number of times to perform the halo exchange for timing
    opts.addOpt< int >( "nexchanges", "Number of ghost-halo exchange iterations to perform. Default=10",
                        &num_max_exchange );
 
    opts.parseCommandLine( argc, argv );
}
 
void RuntimeContext::timer_push( const std::string& operation )
{
    mTimerOps = mTimer.time_since_birth();
    mOpName   = operation;
}
 
void RuntimeContext::timer_pop( const int nruns )
{
    double locElapsed = mTimer.time_since_birth() - mTimerOps;
    double avgElapsed = 0;
    double maxElapsed = 0;
    MPI_Reduce( &locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, parallel_communicator->comm() );
    MPI_Reduce( &locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, parallel_communicator->comm() );
    if( proc_id == 0 )
    {
        avgElapsed /= num_procs;
        if( nruns > 1 )
            std::cout << "[LOG] Time taken to " << mOpName.c_str() << ", averaged over " << nruns
                      << " runs : max = " << maxElapsed / nruns << ", avg = " << avgElapsed / nruns << "\n";
        else
            std::cout << "[LOG] Time taken to " << mOpName.c_str() << " : max = " << maxElapsed
                      << ", avg = " << avgElapsed << "\n";
 
        last_counter = maxElapsed / nruns;
    }
    mOpName.clear();
}
 
double RuntimeContext::last_elapsed() const
{
    return last_counter;
}
 
moab::ErrorCode RuntimeContext::create_sv_tags( moab::Tag& tagScalar, moab::Tag& tagVector, moab::Range& entities )
{
    // Get element (centroid) coordinates so that we can evaluate some arbitrary data
    std::vector< double > entCoords = compute_centroids( entities );  // [entities * [lon, lat]]
 
    if( proc_id == 0 ) std::cout << "> Getting scalar tag handle " << scalar_tagname << "..." << std::endl;
    double defSTagValue = -1.0;
    bool createdTScalar = false;
    // Get or create the scalar tag: default name = "scalar_variable"
    // Type: double, Components: 1, Layout: Dense (all entities potentially), Default: -1.0
    runchk( moab_interface.tag_get_handle( scalar_tagname.c_str(), 1, moab::MB_TYPE_DOUBLE, tagScalar,
                                           moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &defSTagValue, &createdTScalar ),
            "Retrieving scalar tag handle failed" );
 
    // we expect to create a new tag -- fail if Tag already exists since we do not want to overwrite data
    assert( createdTScalar );
    // set the data for scalar tag with an analytical Spherical Harmonic function
    {
        std::vector< double > tagValues( entities.size(), -1.0 );
        std::generate( tagValues.begin(), tagValues.end(), [=, &entCoords]() {
            static int index = 0;
            const int offset = index++ * 2;
            return evaluate_function( entCoords[offset], entCoords[offset + 1] );
        } );
        // Set local scalar tag data for exchange
        runchk( moab_interface.tag_set_data( tagScalar, entities, tagValues.data() ),
                "Setting scalar tag data failed" );
    }
 
    if( proc_id == 0 ) std::cout << "> Getting vector tag handle " << vector_tagname << "..." << std::endl;
    std::vector< double > defVTagValue( vector_length, -1.0 );
    bool createdTVector = false;
    // Get or create the scalar tag: default name = "vector_variable"
    // Type: double, Components: vector_length, Layout: Dense (all entities potentially), Default: [-1.0,..]
    runchk( moab_interface.tag_get_handle( vector_tagname.c_str(), vector_length, moab::MB_TYPE_DOUBLE, tagVector,
                                           moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, defVTagValue.data(),
                                           &createdTVector ),
            "Retrieving vector tag handle failed" );
 
    // we expect to create a new tag -- fail if Tag already exists since we do not want to overwrite data
    assert( createdTVector );
    // set the data for vector tag with an analytical Spherical Harmonic function
    // with an optional scaling for each component; just to make it look different :-)
    {
        const int veclength = vector_length;
        std::vector< double > tagValues( entities.size() * veclength, -1.0 );
        std::generate( tagValues.begin(), tagValues.end(), [=, &entCoords]() {
            static int index = 0;
            const int offset = ( index++ / veclength ) * 2;
            return this->evaluate_function( entCoords[offset], entCoords[offset + 1], 2, ( index % veclength + 1.0 ) );
        } );
        // Set local tag data for exchange
        runchk( moab_interface.tag_set_data( tagVector, entities, tagValues.data() ),
                "Setting vector tag data failed" );
    }
 
    return moab::MB_SUCCESS;
}
 
moab::ErrorCode RuntimeContext::load_file( bool load_ghosts )
{
    /// Parallel Read options:
    ///   PARALLEL = type {READ_PART} : Read on all tasks
    ///   PARTITION_METHOD = RCBZOLTAN : Use Zoltan partitioner to compute an online partition and redistribute on the
    ///   fly PARTITION = PARALLEL_PARTITION : Partition as you read based on part information stored in h5m file
    ///   PARALLEL_RESOLVE_SHARED_ENTS : Communicate to all processors to get the shared adjacencies
    ///   consistently in parallel
    ///   PARALLEL_GHOSTS : a.b.c
    ///                   : a = 2 - highest dimension of entities (2D in this case)
    ///                   : b = 1 - dimension of entities to calculate adjacencies (vertex=0, edges=1)
    ///                   : c = 3 - number of ghost layers needed (3 in this case)
    std::string read_options   = "DEBUG_IO=0;";
    std::string::size_type idx = input_filename.rfind( '.' );
    if( num_procs > 1 && idx != std::string::npos )
    {
        std::string extension = input_filename.substr( idx + 1 );
        if( !extension.compare( "nc" ) )
            // PARTITION_METHOD= [RCBZOLTAN, TRIVIAL]
            read_options += "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
                            "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
        else if( !extension.compare( "h5m" ) )
            read_options +=
                "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;"
                "PARALLEL_RESOLVE_SHARED_ENTS;" +
                ( load_ghosts ? "PARALLEL_THIN_GHOST_LAYER;PARALLEL_GHOSTS=2.1." + std::to_string( ghost_layers ) + ";"
                              : "" );
        else
        {
            std::cout << "Error unsupported file type (only h5m and nc) for this example: " << input_filename
                      << std::endl;
            return moab::MB_UNSUPPORTED_OPERATION;
        }
    }
 
    // Load the file from disk with given read options in parallel and associate all entities to fileset
    return moab_interface.load_file( input_filename.c_str(), &fileset, read_options.c_str() );
}
 
std::vector< double > RuntimeContext::compute_centroids( const moab::Range& entities ) const
{
    double node[3];
    std::vector< double > eCentroids( entities.size() * 2 );  // [lon, lat]
    size_t offset = 0;
    for( auto entity : entities )
    {
        // Get the element coordinates (centroid) on the real mesh
        runchk_cont( moab_interface.get_coords( &entity, 1, node ), "Getting entity coordinates failed" );
 
        // scale by magnitude so that element is on unit sphere
        double magnitude = std::sqrt( node[0] * node[0] + node[1] * node[1] + node[2] * node[2] );
        node[0] /= magnitude;
        node[1] /= magnitude;
        node[2] /= magnitude;
 
        // compute the spherical transformation onto unit sphere
        eCentroids[offset] = atan2( node[1], node[0] );
        if( eCentroids[offset] < 0.0 ) eCentroids[offset] += 2.0 * M_PI;
        eCentroids[offset + 1] = asin( node[2] );
 
        offset += 2;  // increment the offset
    }
    // return centroid list for elements
    return eCentroids;
}