Cppcheck report - MOAB: src/ScdInterface.cpp

#include "moab/ScdInterface.hpp"
#include "moab/Core.hpp"
#include "SequenceManager.hpp"
#include "EntitySequence.hpp"
#include "StructuredElementSeq.hpp"
#include "VertexSequence.hpp"
#include "ScdVertexData.hpp"
#include "MBTagConventions.hpp"
#ifdef MOAB_HAVE_MPI
#include "moab/ParallelComm.hpp"
#include "moab/TupleList.hpp"
#include "moab/gs.hpp"
#endif
#include <cassert>
#include <iostream>
#include <functional>

#define ERRORR( rval, str )          \
    {                                \
        if( MB_SUCCESS != ( rval ) ) \
        {                            \
            std::cerr << ( str );    \
            return rval;             \
        }                            \
    }

namespace moab
{

const char* ScdParData::PartitionMethodNames[] = { "alljorkori", "alljkbal", "sqij",      "sqjk",
                                                   "sqijk",      "trivial",  "rcbzoltan", "nopart" };

ScdInterface::ScdInterface( Interface* imp, bool boxes )
    : mbImpl( imp ), searchedBoxes( false ), boxPeriodicTag( 0 ), boxDimsTag( 0 ), globalBoxDimsTag( 0 ),
      partMethodTag( 0 ), boxSetTag( 0 )
{
    if( boxes ) find_boxes( scdBoxes );
}

// Destructor
ScdInterface::~ScdInterface()
{
    std::vector< ScdBox* > tmp_boxes;
    tmp_boxes.swap( scdBoxes );

    for( std::vector< ScdBox* >::iterator rit = tmp_boxes.begin(); rit != tmp_boxes.end(); ++rit )
        delete *rit;

    if( box_set_tag( false ) ) mbImpl->tag_delete( box_set_tag() );
}

Interface* ScdInterface::impl() const
{
    return mbImpl;
}

ErrorCode ScdInterface::find_boxes( std::vector< ScdBox* >& scd_boxes )
{
    Range tmp_boxes;
    ErrorCode rval = find_boxes( tmp_boxes );
    if( MB_SUCCESS != rval ) return rval;

    for( Range::iterator rit = tmp_boxes.begin(); rit != tmp_boxes.end(); ++rit )
    {
        ScdBox* tmp_box = get_scd_box( *rit );
        if( tmp_box )
            scd_boxes.push_back( tmp_box );
        else
            rval = MB_FAILURE;
    }

    return rval;
}

ErrorCode ScdInterface::find_boxes( Range& scd_boxes )
{
    ErrorCode rval = MB_SUCCESS;
    box_dims_tag();
    Range boxes;
    if( !searchedBoxes )
    {
        rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &boxDimsTag, NULL, 1, boxes, Interface::UNION );
        searchedBoxes = true;
        if( !boxes.empty() )
        {
            scdBoxes.resize( boxes.size() );
            rval        = mbImpl->tag_get_data( boxSetTag, boxes, &scdBoxes[0] );
            ScdBox* dum = NULL;
            // std::remove_if(scdBoxes.begin(), scdBoxes.end(),
            // std::bind2nd(std::equal_to<ScdBox*>(), dum) ) ;
            std::remove_if( scdBoxes.begin(), scdBoxes.end(),<--- Return value of function std::remove_if() is not used.
<--- Return value of std::remove_if() ignored. Elements remain in container. [+]
The return value of std::remove_if() is ignored. This function returns an iterator to the end of the range containing those elements that should be kept. Elements past new end remain valid but with unspecified values. Use the erase method of the container to delete them.
                            std::bind( std::equal_to< ScdBox* >(), std::placeholders::_1, dum ) );
            // https://stackoverflow.com/questions/32739018/a-replacement-for-stdbind2nd
        }
    }

    for( std::vector< ScdBox* >::iterator vit = scdBoxes.begin(); vit != scdBoxes.end(); ++vit )
        scd_boxes.insert( ( *vit )->box_set() );

    return rval;
}

ScdBox* ScdInterface::get_scd_box( EntityHandle eh )
{
    ScdBox* scd_box = NULL;
    if( !box_set_tag( false ) ) return scd_box;

    mbImpl->tag_get_data( box_set_tag(), &eh, 1, &scd_box );
    return scd_box;
}

ErrorCode ScdInterface::construct_box( HomCoord low,
                                       HomCoord high,
                                       const double* const coords,
                                       unsigned int num_coords,
                                       ScdBox*& new_box,
                                       int* const lperiodic,
                                       ScdParData* par_data,
                                       bool assign_gids,
                                       int tag_shared_ents )
{
    // create a rectangular structured mesh block
    ErrorCode rval;

    int tmp_lper[3] = { 0, 0, 0 };
    if( lperiodic ) std::copy( lperiodic, lperiodic + 3, tmp_lper );

#ifndef MOAB_HAVE_MPI
    if( -1 != tag_shared_ents ) ERRORR( MB_FAILURE, "Parallel capability requested but MOAB not compiled parallel." );
    if( -1 == tag_shared_ents && !assign_gids ) assign_gids = true;  // need to assign gids in order to tag shared verts
#else
    if( par_data && low == high && ScdParData::NOPART != par_data->partMethod )
    {
        // requesting creation of parallel mesh, so need to compute partition
        if( !par_data->pComm )
        {
            // this is a really boneheaded way to have to create a PC
            par_data->pComm = ParallelComm::get_pcomm( mbImpl, 0 );
            if( NULL == par_data->pComm ) par_data->pComm = new ParallelComm( mbImpl, MPI_COMM_WORLD );
        }
        int ldims[6];
        rval = compute_partition( par_data->pComm->size(), par_data->pComm->rank(), *par_data, ldims, tmp_lper,
                                  par_data->pDims );ERRORR( rval, "Error returned from compute_partition." );
        low.set( ldims[0], ldims[1], ldims[2] );
        high.set( ldims[3], ldims[4], ldims[5] );
        if( par_data->pComm->get_debug_verbosity() > 0 )
        {
            std::cout << "Proc " << par_data->pComm->rank() << ": " << *par_data;
            std::cout << "Proc " << par_data->pComm->rank() << " local dims: " << low << "-" << high << std::endl;
        }
    }
#endif

    HomCoord tmp_size = high - low + HomCoord( 1, 1, 1, 0 );
    if( ( tmp_size[1] && num_coords && (int)num_coords < tmp_size[0] ) ||
        ( tmp_size[2] && num_coords && (int)num_coords < tmp_size[0] * tmp_size[1] ) )
        return MB_FAILURE;

    rval = create_scd_sequence( low, high, MBVERTEX, 0, new_box );ERRORR( rval, "Trouble creating scd vertex sequence." );

    // set the vertex coordinates
    double *xc, *yc, *zc;
    rval = new_box->get_coordinate_arrays( xc, yc, zc );ERRORR( rval, "Couldn't get vertex coordinate arrays." );

    if( coords && num_coords )
    {
        unsigned int i = 0;
        for( int kl = low[2]; kl <= high[2]; kl++ )
        {
            for( int jl = low[1]; jl <= high[1]; jl++ )
            {
                for( int il = low[0]; il <= high[0]; il++ )
                {
                    xc[i] = coords[3 * i];
                    if( new_box->box_size()[1] ) yc[i] = coords[3 * i + 1];
                    if( new_box->box_size()[2] ) zc[i] = coords[3 * i + 2];
                    i++;
                }
            }
        }
    }
    else
    {
        unsigned int i = 0;
        for( int kl = low[2]; kl <= high[2]; kl++ )
        {
            for( int jl = low[1]; jl <= high[1]; jl++ )
            {
                for( int il = low[0]; il <= high[0]; il++ )
                {
                    xc[i] = (double)il;
                    if( new_box->box_size()[1] )
                        yc[i] = (double)jl;
                    else
                        yc[i] = 0.0;
                    if( new_box->box_size()[2] )
                        zc[i] = (double)kl;
                    else
                        zc[i] = 0.0;
                    i++;
                }
            }
        }
    }

    // create element sequence
    Core* mbcore             = dynamic_cast< Core* >( mbImpl );
    SequenceManager* seq_mgr = mbcore->sequence_manager();

    EntitySequence* tmp_seq;
    EntityHandle start_ent;

    // construct the sequence
    EntityType this_tp = MBHEX;
    if( 1 >= tmp_size[2] ) this_tp = MBQUAD;
    if( 1 >= tmp_size[2] && 1 >= tmp_size[1] ) this_tp = MBEDGE;
    rval = seq_mgr->create_scd_sequence( low, high, this_tp, 0, start_ent, tmp_seq, tmp_lper );ERRORR( rval, "Trouble creating scd element sequence." );

    new_box->elem_seq( tmp_seq );
    new_box->start_element( start_ent );

    // add vertex seq to element seq, forward orientation, unity transform
    rval = new_box->add_vbox( new_box,
                              // p1: imin,jmin
                              low, low,
                              // p2: imax,jmin
                              low + HomCoord( 1, 0, 0 ), low + HomCoord( 1, 0, 0 ),
                              // p3: imin,jmax
                              low + HomCoord( 0, 1, 0 ), low + HomCoord( 0, 1, 0 ) );ERRORR( rval, "Error constructing structured element sequence." );

    // add the new hexes to the scd box set; vertices were added in call to create_scd_sequence
    Range tmp_range( new_box->start_element(), new_box->start_element() + new_box->num_elements() - 1 );
    rval = mbImpl->add_entities( new_box->box_set(), tmp_range );ERRORR( rval, "Couldn't add new hexes to box set." );

    if( par_data ) new_box->par_data( *par_data );

    if( assign_gids )
    {
        rval = assign_global_ids( new_box );ERRORR( rval, "Trouble assigning global ids" );
    }

#ifdef MOAB_HAVE_MPI
    if( par_data && -1 != tag_shared_ents )
    {
        rval = tag_shared_vertices( par_data->pComm, new_box );
    }
#endif

    return MB_SUCCESS;
}

ErrorCode ScdInterface::assign_global_ids( ScdBox* box )
{
    // Get a ptr to global id memory
    void* data;
    int count   = 0;
    Tag gid_tag = mbImpl->globalId_tag();
    Range tmp_range( box->start_vertex(), box->start_vertex() + box->num_vertices() );
    ErrorCode rval = mbImpl->tag_iterate( gid_tag, tmp_range.begin(), tmp_range.end(), count, data );ERRORR( rval, "Failed to get tag iterator." );
    assert( count == box->num_vertices() );
    int* gid_data = (int*)data;
    int di        = box->par_data().gDims[3] - box->par_data().gDims[0] + 1;
    int dj        = box->par_data().gDims[4] - box->par_data().gDims[1] + 1;

    for( int kl = box->box_dims()[2]; kl <= box->box_dims()[5]; kl++ )
    {
        for( int jl = box->box_dims()[1]; jl <= box->box_dims()[4]; jl++ )
        {
            for( int il = box->box_dims()[0]; il <= box->box_dims()[3]; il++ )
            {
                int itmp =
                    ( !box->locally_periodic()[0] && box->par_data().gPeriodic[0] && il == box->par_data().gDims[3]
                          ? box->par_data().gDims[0]
                          : il );
                *gid_data = ( -1 != kl ? kl * di * dj : 0 ) + jl * di + itmp + 1;
                gid_data++;
            }
        }
    }

    return MB_SUCCESS;
}

ErrorCode ScdInterface::create_scd_sequence( const HomCoord& low,
                                             const HomCoord& high,
                                             EntityType tp,
                                             int starting_id,
                                             ScdBox*& new_box,
                                             int* is_periodic )
{
    HomCoord tmp_size = high - low + HomCoord( 1, 1, 1, 0 );
    if( ( tp == MBHEX && 1 >= tmp_size[2] ) || ( tp == MBQUAD && 1 >= tmp_size[1] ) ||
        ( tp == MBEDGE && 1 >= tmp_size[0] ) )
        return MB_TYPE_OUT_OF_RANGE;

    Core* mbcore = dynamic_cast< Core* >( mbImpl );
    assert( mbcore != NULL );
    SequenceManager* seq_mgr = mbcore->sequence_manager();

    EntitySequence* tmp_seq;
    EntityHandle start_ent, scd_set;

    // construct the sequence
    ErrorCode rval = seq_mgr->create_scd_sequence( low, high, tp, starting_id, start_ent, tmp_seq, is_periodic );
    if( MB_SUCCESS != rval ) return rval;

    // create the set for this rectangle
    rval = create_box_set( low, high, scd_set );
    if( MB_SUCCESS != rval ) return rval;

    // make the ScdBox
    new_box = new ScdBox( this, scd_set, tmp_seq );
    if( !new_box ) return MB_FAILURE;

    // set the start vertex/element
    Range new_range;
    if( MBVERTEX == tp )
    {
        new_range.insert( start_ent, start_ent + new_box->num_vertices() - 1 );
    }
    else
    {
        new_range.insert( start_ent, start_ent + new_box->num_elements() - 1 );
    }

    // put the entities in the box set
    rval = mbImpl->add_entities( scd_set, new_range );
    if( MB_SUCCESS != rval ) return rval;

    // tag the set with the box
    rval = mbImpl->tag_set_data( box_set_tag(), &scd_set, 1, &new_box );
    if( MB_SUCCESS != rval ) return rval;

    return MB_SUCCESS;
}

ErrorCode ScdInterface::create_box_set( const HomCoord& low,
                                        const HomCoord& high,
                                        EntityHandle& scd_set,
                                        int* is_periodic )
{
    // create the set and put the entities in it
    ErrorCode rval = mbImpl->create_meshset( MESHSET_SET, scd_set );
    if( MB_SUCCESS != rval ) return rval;

    // tag the set with parameter extents
    int boxdims[6];
    for( int i = 0; i < 3; i++ )
        boxdims[i] = low[i];
    for( int i = 0; i < 3; i++ )
        boxdims[3 + i] = high[i];
    rval = mbImpl->tag_set_data( box_dims_tag(), &scd_set, 1, boxdims );
    if( MB_SUCCESS != rval ) return rval;

    if( is_periodic )
    {
        rval = mbImpl->tag_set_data( box_periodic_tag(), &scd_set, 1, is_periodic );
        if( MB_SUCCESS != rval ) return rval;
    }

    return rval;
}

Tag ScdInterface::box_periodic_tag( bool create_if_missing )
{
    // Reset boxPeriodicTag in case it has been deleted (e.g. by Core::clean_up_failed_read)
    if( boxPeriodicTag )
    {
        std::string tag_name;
        if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxPeriodicTag, tag_name ) ) boxPeriodicTag = NULL;
    }

    if( boxPeriodicTag || !create_if_missing ) return boxPeriodicTag;

    ErrorCode rval =
        mbImpl->tag_get_handle( "BOX_PERIODIC", 3, MB_TYPE_INTEGER, boxPeriodicTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != rval ) return 0;
    return boxPeriodicTag;
}

Tag ScdInterface::box_dims_tag( bool create_if_missing )
{
    // Reset boxDimsTag in case it has been deleted (e.g. by clean_up_failed_read)
    if( boxDimsTag )
    {
        std::string tag_name;
        if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxDimsTag, tag_name ) ) boxDimsTag = NULL;
    }

    if( boxDimsTag || !create_if_missing ) return boxDimsTag;

    ErrorCode rval = mbImpl->tag_get_handle( "BOX_DIMS", 6, MB_TYPE_INTEGER, boxDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != rval ) return 0;
    return boxDimsTag;
}

Tag ScdInterface::global_box_dims_tag( bool create_if_missing )<--- The function 'global_box_dims_tag' is never used.
{
    // Reset globalBoxDimsTag in case it has been deleted (e.g. by Core::clean_up_failed_read)
    if( globalBoxDimsTag )
    {
        std::string tag_name;
        if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( globalBoxDimsTag, tag_name ) ) globalBoxDimsTag = NULL;
    }

    if( globalBoxDimsTag || !create_if_missing ) return globalBoxDimsTag;

    ErrorCode rval =
        mbImpl->tag_get_handle( "GLOBAL_BOX_DIMS", 6, MB_TYPE_INTEGER, globalBoxDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != rval ) return 0;
    return globalBoxDimsTag;
}

Tag ScdInterface::part_method_tag( bool create_if_missing )
{
    // Reset partMethodTag in case it has been deleted (e.g. by Core::clean_up_failed_read)
    if( partMethodTag )
    {
        std::string tag_name;
        if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( partMethodTag, tag_name ) ) partMethodTag = NULL;
    }

    if( partMethodTag || !create_if_missing ) return partMethodTag;

    ErrorCode rval =
        mbImpl->tag_get_handle( "PARTITION_METHOD", 1, MB_TYPE_INTEGER, partMethodTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != rval ) return 0;
    return partMethodTag;
}

Tag ScdInterface::box_set_tag( bool create_if_missing )
{
    // Reset boxSetTag in case it has been deleted (e.g. by Core::clean_up_failed_read)
    if( boxSetTag )
    {
        std::string tag_name;
        if( MB_TAG_NOT_FOUND == mbImpl->tag_get_name( boxSetTag, tag_name ) ) boxSetTag = NULL;
    }

    if( boxSetTag || !create_if_missing ) return boxSetTag;

    ErrorCode rval = mbImpl->tag_get_handle( "__BOX_SET", sizeof( ScdBox* ), MB_TYPE_OPAQUE, boxSetTag,
                                             MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != rval ) return 0;
    return boxSetTag;
}

//! Remove the box from the list on ScdInterface
ErrorCode ScdInterface::remove_box( ScdBox* box )
{
    std::vector< ScdBox* >::iterator vit = std::find( scdBoxes.begin(), scdBoxes.end(), box );
    if( vit != scdBoxes.end() )
    {
        scdBoxes.erase( vit );
        return MB_SUCCESS;
    }
    else
        return MB_FAILURE;
}

//! Add the box to the list on ScdInterface
ErrorCode ScdInterface::add_box( ScdBox* box )
{
    scdBoxes.push_back( box );
    return MB_SUCCESS;
}

ErrorCode ScdInterface::get_boxes( std::vector< ScdBox* >& boxes )<--- The function 'get_boxes' is never used.
{
    std::copy( scdBoxes.begin(), scdBoxes.end(), std::back_inserter( boxes ) );
    return MB_SUCCESS;
}

ScdBox::ScdBox( ScdInterface* impl, EntityHandle bset, EntitySequence* seq1, EntitySequence* seq2 )
    : scImpl( impl ), boxSet( bset ), vertDat( NULL ), elemSeq( NULL ), startVertex( 0 ), startElem( 0 )
{
    for( int i = 0; i < 6; i++ )
        boxDims[i] = 0;
    for( int i = 0; i < 3; i++ )
        locallyPeriodic[i] = false;
    VertexSequence* vseq = dynamic_cast< VertexSequence* >( seq1 );
    if( vseq ) vertDat = dynamic_cast< ScdVertexData* >( vseq->data() );
    if( vertDat )
    {
        // retrieve the parametric space
        for( int i = 0; i < 3; i++ )
        {
            boxDims[i]     = vertDat->min_params()[i];
            boxDims[3 + i] = vertDat->max_params()[i];
        }
        startVertex = vertDat->start_handle();
    }
    else if( impl->boxDimsTag )
    {
        // look for parametric space info on set
        ErrorCode rval = impl->mbImpl->tag_get_data( impl->boxDimsTag, &bset, 1, boxDims );
        if( MB_SUCCESS == rval )
        {
            Range verts;
            impl->mbImpl->get_entities_by_dimension( bset, 0, verts );
            if( !verts.empty() ) startVertex = *verts.begin();
        }
    }

    elemSeq = dynamic_cast< StructuredElementSeq* >( seq2 );
    if( !elemSeq ) elemSeq = dynamic_cast< StructuredElementSeq* >( seq1 );

    if( elemSeq )
    {
        if( vertDat )
        {
            // check the parametric space to make sure it's consistent
            assert( elemSeq->sdata()->min_params() == HomCoord( boxDims, 3 ) &&
<--- Assert statement calls a function which may have desired side effects: 'sdata'. [+]
Non-pure function: 'sdata' is called inside assert statement. Assert statements are removed from release builds so the code inside assert statement is not executed. If the code is needed also in release builds, this is a bug.
                    ( elemSeq->sdata()->max_params() + HomCoord( 1, 1, 1 ) ) == HomCoord( boxDims, 3 ) );
<--- Assert statement calls a function which may have desired side effects: 'sdata'. [+]
Non-pure function: 'sdata' is called inside assert statement. Assert statements are removed from release builds so the code inside assert statement is not executed. If the code is needed also in release builds, this is a bug.
        }
        else
        {
            // get the parametric space from the element sequence
            for( int i = 0; i < 3; i++ )
            {
                boxDims[i]     = elemSeq->sdata()->min_params()[i];
                boxDims[3 + i] = elemSeq->sdata()->max_params()[i];
            }
        }

        startElem = elemSeq->start_handle();
    }
    else
    {
        Range elems;
        impl->mbImpl->get_entities_by_dimension(
            bset, ( boxDims[2] == boxDims[5] ? ( boxDims[1] == boxDims[4] ? 1 : 2 ) : 3 ), elems );
        if( !elems.empty() ) startElem = *elems.begin();
        // call the following w/o looking at return value, since it doesn't really need to be there
        if( impl->boxPeriodicTag ) impl->mbImpl->tag_get_data( impl->boxPeriodicTag, &bset, 1, locallyPeriodic );
    }

    assert( vertDat || elemSeq || boxDims[0] != boxDims[3] || boxDims[1] != boxDims[4] || boxDims[2] != boxDims[5] );

    boxSize     = HomCoord( boxDims + 3, 3 ) - HomCoord( boxDims, 3 ) + HomCoord( 1, 1, 1 );
    boxSizeIJ   = ( boxSize[1] ? boxSize[1] : 1 ) * boxSize[0];
    boxSizeIM1  = boxSize[0] - ( locallyPeriodic[0] ? 0 : 1 );
    boxSizeIJM1 = ( boxSize[1] ? ( boxSize[1] - ( locallyPeriodic[1] ? 0 : 1 ) ) : 1 ) * boxSizeIM1;

    scImpl->add_box( this );
}

ScdBox::~ScdBox()
{
    // Reset the tag on the set
    if( boxSet )
    {
        // It is possible that the box set entity has been deleted (e.g. by
        // Core::clean_up_failed_read)
        Core* mbcore = dynamic_cast< Core* >( scImpl->mbImpl );
        assert( mbcore != NULL );
        if( mbcore->is_valid( boxSet ) )
        {
            ScdBox* tmp_ptr = NULL;
            scImpl->mbImpl->tag_set_data( scImpl->box_set_tag(), &boxSet, 1, &tmp_ptr );
        }
        else
            boxSet = 0;
    }

    scImpl->remove_box( this );
}

EntityHandle ScdBox::get_vertex_from_seq( int i, int j, int k ) const
{
    assert( elemSeq );
    return elemSeq->get_vertex( i, j, k );
}

int ScdBox::box_dimension() const
{
    return ( startElem ? scImpl->mbImpl->dimension_from_handle( startElem ) : -1 );
}

ErrorCode ScdBox::add_vbox( ScdBox* vbox,
                            HomCoord from1,
                            HomCoord to1,
                            HomCoord from2,
                            HomCoord to2,
                            HomCoord from3,
                            HomCoord to3,
                            bool bb_input,
                            const HomCoord& bb_min,
                            const HomCoord& bb_max )
{
    if( !vbox->vertDat ) return MB_FAILURE;
    ScdVertexData* dum_data = dynamic_cast< ScdVertexData* >( vbox->vertDat );
    ErrorCode rval =
        elemSeq->sdata()->add_vsequence( dum_data, from1, to1, from2, to2, from3, to3, bb_input, bb_min, bb_max );
    return rval;
}

bool ScdBox::boundary_complete() const
{
    return elemSeq->boundary_complete();
}

ErrorCode ScdBox::get_coordinate_arrays( double*& xc, double*& yc, double*& zc )
{
    if( !vertDat ) return MB_FAILURE;

    xc = reinterpret_cast< double* >( vertDat->get_sequence_data( 0 ) );
    yc = reinterpret_cast< double* >( vertDat->get_sequence_data( 1 ) );
    zc = reinterpret_cast< double* >( vertDat->get_sequence_data( 2 ) );
    return MB_SUCCESS;
}

ErrorCode ScdBox::get_coordinate_arrays( const double*& xc, const double*& yc, const double*& zc ) const
{
    if( !vertDat ) return MB_FAILURE;
    xc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 0 ) );
    yc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 1 ) );
    zc = reinterpret_cast< const double* >( vertDat->get_sequence_data( 2 ) );
    return MB_SUCCESS;
}

ErrorCode ScdBox::vert_dat( ScdVertexData* vert_dt )<--- The function 'vert_dat' is never used.
{
    vertDat = vert_dt;
    return MB_SUCCESS;
}

ErrorCode ScdBox::elem_seq( EntitySequence* elem_sq )
{
    elemSeq = dynamic_cast< StructuredElementSeq* >( elem_sq );
    if( elemSeq ) elemSeq->is_periodic( locallyPeriodic );

    if( locallyPeriodic[0] ) boxSizeIM1 = boxSize[0] - ( locallyPeriodic[0] ? 0 : 1 );
    if( locallyPeriodic[0] || locallyPeriodic[1] )
        boxSizeIJM1 = ( boxSize[1] ? ( boxSize[1] - ( locallyPeriodic[1] ? 0 : 1 ) ) : 1 ) * boxSizeIM1;

    return ( elemSeq ? MB_SUCCESS : MB_FAILURE );
}

ErrorCode ScdBox::get_params( EntityHandle ent, HomCoord& ijkd ) const
{
    // check first whether this is an intermediate entity, so we know what to do
    int dimension = box_dimension();
    int this_dim  = scImpl->impl()->dimension_from_handle( ent );

    if( ( 0 == this_dim && !vertDat ) || ( this_dim && this_dim == dimension ) )
    {
        assert( elemSeq );
        return elemSeq->get_params( ent, ijkd[0], ijkd[1], ijkd[2] );
    }

    else if( !this_dim && vertDat )
        return vertDat->get_params( ent, ijkd[0], ijkd[1], ijkd[2] );

    else
        return MB_NOT_IMPLEMENTED;
}

//! Get the entity of specified dimension adjacent to parametric element
/**
 * \param dim Dimension of adjacent entity being requested
 * \param i Parametric coordinates of cell being evaluated
 * \param j Parametric coordinates of cell being evaluated
 * \param k Parametric coordinates of cell being evaluated
 * \param dir Direction (0, 1, or 2), for getting adjacent edges (2d, 3d) or faces (3d)
 * \param ent EntityHandle of adjacent entity
 * \param create_if_missing If true, creates the entity if it doesn't already exist
 */
ErrorCode ScdBox::get_adj_edge_or_face( int dim,
                                        int i,
                                        int j,
                                        int k,
                                        int dir,
                                        EntityHandle& ent,
                                        bool create_if_missing ) const
{
    // describe connectivity of sub-element in static array
    // subconnect[dim-1][dir][numv][ijk] where dimensions are:
    // [dim-1]: dim=1 or 2, so this is 0 or 1
    // [dir]: one of 0..2, for ijk directions in a hex
    // [numv]: number of vertices describing sub entity = 2*dim <= 4
    // [ijk]: 3 values for i, j, k
    int subconnect[2][3][4][3] = { { { { 0, 0, 0 }, { 1, 0, 0 }, { -1, -1, -1 }, { -1, -1, -1 } },    // i edge
                                     { { 0, 0, 0 }, { 0, 1, 0 }, { -1, -1, -1 }, { -1, -1, -1 } },    // j edge
                                     { { 0, 0, 0 }, { 0, 0, 1 }, { -1, -1, -1 }, { -1, -1, -1 } } },  // k edge

                                   { { { 0, 0, 0 }, { 0, 1, 0 }, { 0, 1, 1 }, { 0, 0, 1 } },      // i face
                                     { { 0, 0, 0 }, { 1, 0, 0 }, { 1, 0, 1 }, { 0, 0, 1 } },      // j face
                                     { { 0, 0, 0 }, { 1, 0, 0 }, { 1, 1, 0 }, { 0, 1, 0 } } } };  // k face

    // check proper input dimensions and lower bound
    if( dim < 1 || dim > 2 || i < boxDims[0] || j < boxDims[1] || k < boxDims[2] ) return MB_FAILURE;

    // now check upper bound; parameters must be <= upper corner, since edges/faces
    // follow element parameterization, not vertex parameterization
    else if( ( boxDims[3] != boxDims[0] && i > ( locallyPeriodic[0] ? boxDims[3] + 1 : boxDims[3] ) ) ||
             ( boxDims[4] != boxDims[1] && j > ( locallyPeriodic[1] ? boxDims[4] + 1 : boxDims[4] ) ) ||
             ( boxDims[5] != boxDims[2] && k > boxDims[5] ) )
        return MB_FAILURE;

    // get the vertices making up this entity
    EntityHandle verts[4];
    for( int ind = 0; ind < 2 * dim; ind++ )
    {
        int i1 = i + subconnect[dim - 1][dir][ind][0];
        int j1 = j + subconnect[dim - 1][dir][ind][1];
        // if periodic in i and i1 is boxDims[3]+1, wrap around
        if( locallyPeriodic[0] && i1 == boxDims[3] + 1 ) i1 = boxDims[0];
        // if periodic in j and j1 is boxDims[4]+1, wrap around
        if( locallyPeriodic[1] && j1 == boxDims[4] + 1 ) j1 = boxDims[1];
        verts[ind] = get_vertex( i1, j1, k + subconnect[dim - 1][dir][ind][2] );
        if( !verts[ind] ) return MB_FAILURE;
    }

    Range ents;
    ErrorCode rval = scImpl->impl()->get_adjacencies( verts, 2 * dim, dim, false, ents );
    if( MB_SUCCESS != rval ) return rval;

    if( ents.size() > 1 )
        return MB_FAILURE;

    else if( ents.size() == 1 )
    {
        ent = *ents.begin();
    }
    else if( create_if_missing )
        rval = scImpl->impl()->create_element( ( 1 == dim ? MBEDGE : MBQUAD ), verts, 2 * dim, ent );

    return rval;
}

#ifndef MOAB_HAVE_MPI
ErrorCode ScdInterface::tag_shared_vertices( ParallelComm*, ScdBox* )
{
    return MB_FAILURE;
#else
ErrorCode ScdInterface::tag_shared_vertices( ParallelComm* pcomm, ScdBox* box )
{
    EntityHandle seth = box->box_set();

    // check the # ents in the box against the num in the set, to make sure it's only 1 box;
    // reuse tmp_range
    Range tmp_range;
    ErrorCode rval = mbImpl->get_entities_by_dimension( seth, box->box_dimension(), tmp_range );
    if( MB_SUCCESS != rval ) return rval;
    if( box->num_elements() != (int)tmp_range.size() ) return MB_FAILURE;

    const int* gdims = box->par_data().gDims;
    if( ( gdims[0] == gdims[3] && gdims[1] == gdims[4] && gdims[2] == gdims[5] ) || -1 == box->par_data().partMethod )
        return MB_FAILURE;

    // ok, we have a partitioned box; get the vertices shared with other processors
    std::vector< int > procs, offsets, shared_indices;
    rval = get_shared_vertices( pcomm, box, procs, offsets, shared_indices );
    if( MB_SUCCESS != rval ) return rval;

    // post receives for start handles once we know how many to look for
    std::vector< MPI_Request > recv_reqs( procs.size(), MPI_REQUEST_NULL ), send_reqs( procs.size(), MPI_REQUEST_NULL );
    std::vector< EntityHandle > rhandles( 4 * procs.size() ), shandles( 4 );
    for( unsigned int i = 0; i < procs.size(); i++ )
    {
        int success = MPI_Irecv( (void*)&rhandles[4 * i], 4 * sizeof( EntityHandle ), MPI_UNSIGNED_CHAR, procs[i], 1,
                                 pcomm->proc_config().proc_comm(), &recv_reqs[i] );
        if( success != MPI_SUCCESS ) return MB_FAILURE;
    }

    // send our own start handles
    shandles[0] = box->start_vertex();
    shandles[1] = 0;
    if( box->box_dimension() == 1 )
    {
        shandles[1] = box->start_element();
        shandles[2] = 0;
        shandles[3] = 0;
    }
    else if( box->box_dimension() == 2 )
    {
        shandles[2] = box->start_element();
        shandles[3] = 0;
    }
    else
    {
        shandles[2] = 0;
        shandles[3] = box->start_element();
    }
    for( unsigned int i = 0; i < procs.size(); i++ )
    {
        int success = MPI_Isend( (void*)&shandles[0], 4 * sizeof( EntityHandle ), MPI_UNSIGNED_CHAR, procs[i], 1,
                                 pcomm->proc_config().proc_comm(), &send_reqs[i] );
        if( success != MPI_SUCCESS ) return MB_FAILURE;
    }

    // receive start handles and save info to a tuple list
    int incoming = procs.size();
    int p, j, k;
    MPI_Status status;
    TupleList shared_data;
    shared_data.initialize( 1, 0, 2, 0, shared_indices.size() / 2 );
    shared_data.enableWriteAccess();

    j = 0;
    k = 0;
    while( incoming )
    {
        int success = MPI_Waitany( procs.size(), &recv_reqs[0], &p, &status );
        if( MPI_SUCCESS != success ) return MB_FAILURE;
        unsigned int num_indices = ( offsets[p + 1] - offsets[p] ) / 2;
        int *lh = &shared_indices[offsets[p]], *rh = lh + num_indices;
        for( unsigned int i = 0; i < num_indices; i++ )
        {
            shared_data.vi_wr[j++] = procs[p];
            shared_data.vul_wr[k++] = shandles[0] + lh[i];
            shared_data.vul_wr[k++] = rhandles[4 * p] + rh[i];
            shared_data.inc_n();
        }
        incoming--;
    }

    // still need to wait for the send requests
    std::vector< MPI_Status > mult_status( procs.size() );
    int success = MPI_Waitall( procs.size(), &send_reqs[0], &mult_status[0] );
    if( MPI_SUCCESS != success )
    {
        MB_SET_ERR( MB_FAILURE, "Failed in waitall in ScdInterface::tag_shared_vertices" );
    }
    // sort by local handle
    TupleList::buffer sort_buffer;
    sort_buffer.buffer_init( shared_indices.size() / 2 );
    shared_data.sort( 1, &sort_buffer );
    sort_buffer.reset();

    // process into sharing data
    std::map< std::vector< int >, std::vector< EntityHandle > > proc_nvecs;
    Range dum;
    rval = pcomm->tag_shared_verts( shared_data, proc_nvecs, dum, 0 );
    if( MB_SUCCESS != rval ) return rval;

    // create interface sets
    rval = pcomm->create_interface_sets( proc_nvecs );
    if( MB_SUCCESS != rval ) return rval;

    // add the box to the PComm's partitionSets
    pcomm->partition_sets().insert( box->box_set() );

    // make sure buffers are allocated for communicating procs
    for( std::vector< int >::iterator pit = procs.begin(); pit != procs.end(); ++pit )
        pcomm->get_buffers( *pit );

    if( pcomm->get_debug_verbosity() > 1 ) pcomm->list_entities( NULL, 1 );

#ifndef NDEBUG
    rval = pcomm->check_all_shared_handles();
    if( MB_SUCCESS != rval ) return rval;
#endif

    return MB_SUCCESS;

#endif
}

ErrorCode ScdInterface::get_neighbor_alljkbal( int np,
                                               int pfrom,
                                               const int* const gdims,
                                               const int* const gperiodic,
                                               const int* const dijk,
                                               int& pto,
                                               int* rdims,
                                               int* facedims,
                                               int* across_bdy )
{
    if( dijk[0] != 0 )
    {
        pto = -1;
        return MB_SUCCESS;
    }

    pto           = -1;
    across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;

    int ldims[6], pijk[3], lperiodic[3];
    ErrorCode rval = compute_partition_alljkbal( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
    if( MB_SUCCESS != rval ) return rval;
    assert( pijk[1] * pijk[2] == np );
    pto        = -1;
    bool bot_j = pfrom< pijk[2], top_j = pfrom > np - pijk[2];
    if( ( 1 == pijk[2] && dijk[2] ) ||                                // 1d in j means no neighbors with dk != 0
        ( !( pfrom % pijk[2] ) && -1 == dijk[2] ) ||                  // at -k bdy
        ( pfrom % pijk[2] == pijk[2] - 1 && 1 == dijk[2] ) ||         // at +k bdy
        ( pfrom < pijk[2] && -1 == dijk[1] && !gperiodic[1] ) ||      // down and not periodic
        ( pfrom >= np - pijk[2] && 1 == dijk[1] && !gperiodic[1] ) )  // up and not periodic
        return MB_SUCCESS;

    pto = pfrom;
    std::copy( ldims, ldims + 6, rdims );
    std::copy( ldims, ldims + 6, facedims );

    if( 0 != dijk[1] )
    {
        pto = ( pto + dijk[1] * pijk[2] + np ) % np;
        assert( pto >= 0 && pto < np );
        int dj = ( gdims[4] - gdims[1] ) / pijk[1], extra = ( gdims[4] - gdims[1] ) % pijk[1];
        if( -1 == dijk[1] )
        {
            facedims[4] = facedims[1];
            if( bot_j )
            {
                // going across periodic lower bdy in j
                rdims[4]      = gdims[4];
                across_bdy[1] = -1;
            }
            else
            {
                rdims[4] = ldims[1];
            }
            rdims[1] = rdims[4] - dj;
            if( pto < extra ) rdims[1]--;
        }
        else
        {
            if( pfrom > np - pijk[2] ) facedims[4] = gdims[1];
            facedims[1] = facedims[4];
            if( top_j )
            {
                // going across periodic upper bdy in j
                rdims[1]      = gdims[1];
                across_bdy[1] = 1;
            }
            else
            {
                rdims[1] = ldims[4];
            }
            rdims[4] = rdims[1] + dj;
            if( pto < extra ) rdims[4]++;
        }
    }
    if( 0 != dijk[2] )
    {
        pto = ( pto + dijk[2] ) % np;
        assert( pto >= 0 && pto < np );
        facedims[2] = facedims[5] = ( -1 == dijk[2] ? facedims[2] : facedims[5] );
        int dk                    = ( gdims[5] - gdims[2] ) / pijk[2];
        if( -1 == dijk[2] )
        {
            facedims[5] = facedims[2];
            rdims[5]    = ldims[2];
            rdims[2]    = rdims[5] - dk;  // never any kextra for alljkbal
        }
        else
        {
            facedims[2] = facedims[5];
            rdims[2]    = ldims[5];
            rdims[5]    = rdims[2] + dk;  // never any kextra for alljkbal
        }
    }

    assert( -1 == pto || ( rdims[0] >= gdims[0] && rdims[3] <= gdims[3] ) );
    assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && bot_j ) ) ) );
    assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) );
    assert( -1 == pto || ( ( facedims[0] >= rdims[0] ||
                             ( gperiodic[0] && rdims[3] == gdims[3] + 1 && facedims[0] == gdims[0] ) ) ) );
    assert( -1 == pto || ( facedims[3] <= rdims[3] ) );
    assert( -1 == pto || ( ( facedims[1] >= rdims[1] ||
                             ( gperiodic[1] && rdims[4] == gdims[4] + 1 && facedims[1] == gdims[1] ) ) ) );
    assert( -1 == pto || ( facedims[4] <= rdims[4] ) );
    assert( -1 == pto || ( facedims[2] >= rdims[2] ) );
    assert( -1 == pto || ( facedims[5] <= rdims[5] ) );
    assert( -1 == pto || ( facedims[0] >= ldims[0] ) );
    assert( -1 == pto || ( facedims[3] <= ldims[3] ) );
    assert( -1 == pto || ( facedims[1] >= ldims[1] ) );
    assert( -1 == pto || ( facedims[4] <= ldims[4] ) );
    assert( -1 == pto || ( facedims[2] >= ldims[2] ) );
    assert( -1 == pto || ( facedims[5] <= ldims[5] ) );

    return MB_SUCCESS;
}

ErrorCode ScdInterface::get_neighbor_sqij( int np,
                                           int pfrom,
                                           const int* const gdims,
                                           const int* const gperiodic,
                                           const int* const dijk,
                                           int& pto,
                                           int* rdims,
                                           int* facedims,
                                           int* across_bdy )
{
    if( dijk[2] != 0 )
    {
        // for sqij, there is no k neighbor, ever
        pto = -1;
        return MB_SUCCESS;
    }

    pto           = -1;
    across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
    int lperiodic[3], pijk[3], ldims[6];
    ErrorCode rval = compute_partition_sqij( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
    if( MB_SUCCESS != rval ) return rval;
    assert( pijk[0] * pijk[1] == np );
    pto        = -1;
    bool top_i = 0, top_j = 0, bot_i = 0, bot_j = 0;
    int ni = pfrom % pijk[0], nj = pfrom / pijk[0];  // row / column number of me
    if( ni == pijk[0] - 1 ) top_i = 1;
    if( nj == pijk[1] - 1 ) top_j = 1;
    if( !ni ) bot_i = 1;
    if( !nj ) bot_j = 1;
    if( ( !gperiodic[0] && bot_i && -1 == dijk[0] ) ||  // left and not periodic
        ( !gperiodic[0] && top_i && 1 == dijk[0] ) ||   // right and not periodic
        ( !gperiodic[1] && bot_j && -1 == dijk[1] ) ||  // bottom and not periodic
        ( !gperiodic[1] && top_j && 1 == dijk[1] ) )    // top and not periodic
        return MB_SUCCESS;

    std::copy( ldims, ldims + 6, facedims );
    std::copy( ldims, ldims + 6, rdims );
    pto   = pfrom;
    int j = gdims[4] - gdims[1], dj = j / pijk[1], jextra = ( gdims[4] - gdims[1] ) % dj, i = gdims[3] - gdims[0],
        di = i / pijk[0], iextra = ( gdims[3] - gdims[0] ) % di;

    if( 0 != dijk[0] )
    {
        pto = ( ni + dijk[0] + pijk[0] ) % pijk[0];  // get pto's ni value
        pto = nj * pijk[0] + pto;                    // then convert to pto
        assert( pto >= 0 && pto < np );
        if( -1 == dijk[0] )
        {
            facedims[3] = facedims[0];
            if( bot_i )
            {
                // going across lower periodic bdy in i
                across_bdy[0] = -1;
                rdims[3]      = gdims[3] + 1;       // +1 because ldims[3] on remote proc is gdims[3]+1
                rdims[0]      = rdims[3] - di - 1;  // -1 to account for rdims[3] being one larger
            }
            else
            {
                rdims[3] = ldims[0];
                rdims[0] = rdims[3] - di;
            }

            if( pto % pijk[0] < iextra ) rdims[0]--;
        }
        else
        {
            if( top_i )
            {
                // going across lower periodic bdy in i
                facedims[3]   = gdims[0];
                across_bdy[0] = 1;
            }
            facedims[0] = facedims[3];
            rdims[0]    = ( top_i ? gdims[0] : ldims[3] );
            rdims[3]    = rdims[0] + di;
            if( pto % pijk[0] < iextra ) rdims[3]++;
            if( gperiodic[0] && ni == pijk[0] - 2 ) rdims[3]++;  // remote proc is top_i and periodic
        }
    }
    if( 0 != dijk[1] )
    {
        pto = ( pto + dijk[1] * pijk[0] + np ) % np;
        assert( pto >= 0 && pto < np );
        if( -1 == dijk[1] )
        {
            facedims[4] = facedims[1];
            if( bot_j )
            {
                // going across lower periodic bdy in j
                rdims[4]      = gdims[4] + 1;       // +1 because ldims[4] on remote proc is gdims[4]+1
                rdims[1]      = rdims[4] - dj - 1;  // -1 to account for gdims[4] being one larger
                across_bdy[1] = -1;
            }
            else
            {
                rdims[4] = ldims[1];
                rdims[1] = rdims[4] - dj;
            }
            if( pto / pijk[0] < jextra ) rdims[1]--;
        }
        else
        {
            if( top_j )
            {
                // going across lower periodic bdy in j
                facedims[4]   = gdims[1];
                rdims[1]      = gdims[1];
                across_bdy[1] = 1;
            }
            else
            {
                rdims[1] = ldims[4];
            }
            facedims[1] = facedims[4];
            rdims[4]    = rdims[1] + dj;
            if( nj + 1 < jextra ) rdims[4]++;
            if( gperiodic[1] && nj == pijk[1] - 2 ) rdims[4]++;  // remote proc is top_j and periodic
        }
    }

    // rdims within gdims
    assert( -1 == pto || ( rdims[0] >= gdims[0] &&
                           ( rdims[3] <= gdims[3] + ( gperiodic[0] && pto % pijk[0] == pijk[0] - 1 ? 1 : 0 ) ) ) );
    assert( -1 == pto || ( rdims[1] >= gdims[1] &&
                           ( rdims[4] <= gdims[4] + ( gperiodic[1] && pto / pijk[0] == pijk[1] - 1 ? 1 : 0 ) ) ) );
    assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) );
    // facedims within rdims
    assert( -1 == pto || ( ( facedims[0] >= rdims[0] ||
                             ( gperiodic[0] && pto % pijk[0] == pijk[0] - 1 && facedims[0] == gdims[0] ) ) ) );
    assert( -1 == pto || ( facedims[3] <= rdims[3] ) );
    assert( -1 == pto || ( ( facedims[1] >= rdims[1] ||
                             ( gperiodic[1] && pto / pijk[0] == pijk[1] - 1 && facedims[1] == gdims[1] ) ) ) );
    assert( -1 == pto || ( facedims[4] <= rdims[4] ) );
    assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) );
    // facedims within ldims
    assert( -1 == pto || ( ( facedims[0] >= ldims[0] || ( top_i && facedims[0] == gdims[0] ) ) ) );
    assert( -1 == pto || ( facedims[3] <= ldims[3] ) );
    assert( -1 == pto || ( ( facedims[1] >= ldims[1] || ( gperiodic[1] && top_j && facedims[1] == gdims[1] ) ) ) );
    assert( -1 == pto || ( facedims[4] <= ldims[4] ) );
    assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) );

    return MB_SUCCESS;
}

ErrorCode ScdInterface::get_neighbor_sqjk( int np,
                                           int pfrom,
                                           const int* const gdims,
                                           const int* const gperiodic,
                                           const int* const dijk,
                                           int& pto,
                                           int* rdims,
                                           int* facedims,
                                           int* across_bdy )
{
    if( dijk[0] != 0 )
    {
        pto = -1;
        return MB_SUCCESS;
    }

    pto           = -1;
    across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
    int pijk[3], lperiodic[3], ldims[6];
    ErrorCode rval = compute_partition_sqjk( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
    if( MB_SUCCESS != rval ) return rval;
    assert( pijk[1] * pijk[2] == np );
    pto        = -1;
    bool top_j = 0, top_k = 0, bot_j = 0, bot_k = 0;
    int nj = pfrom % pijk[1], nk = pfrom / pijk[1];
    if( nj == pijk[1] - 1 ) top_j = 1;
    if( nk == pijk[2] - 1 ) top_k = 1;
    if( !nj ) bot_j = 1;
    if( !nk ) bot_k = 1;
    if( ( !gperiodic[1] && bot_j && -1 == dijk[1] ) ||  // down and not periodic
        ( !gperiodic[1] && top_j && 1 == dijk[1] ) ||   // up and not periodic
        ( bot_k && -1 == dijk[2] ) ||                   // k- bdy
        ( top_k && 1 == dijk[2] ) )                     // k+ bdy
        return MB_SUCCESS;

    std::copy( ldims, ldims + 6, facedims );
    std::copy( ldims, ldims + 6, rdims );
    pto    = pfrom;
    int dj = ( gdims[4] - gdims[1] ) / pijk[1], jextra = ( gdims[4] - gdims[1] ) % dj,
        dk     = ( gdims[5] == gdims[2] ? 0 : ( gdims[5] - gdims[2] ) / pijk[2] ),
        kextra = ( gdims[5] - gdims[2] ) - dk * pijk[2];
    assert( ( dj * pijk[1] + jextra == ( gdims[4] - gdims[1] ) ) &&
            ( dk * pijk[2] + kextra == ( gdims[5] - gdims[2] ) ) );
    if( 0 != dijk[1] )
    {
        pto = ( nj + dijk[1] + pijk[1] ) % pijk[1];  // get pto's ni value
        pto = nk * pijk[1] + pto;                    // then convert to pto
        assert( pto >= 0 && pto < np );
        if( -1 == dijk[1] )
        {
            facedims[4] = facedims[1];
            if( bot_j )
            {
                // going across lower periodic bdy in j
                rdims[4]      = gdims[4] + 1;  // +1 because ldims[4] on remote proc is gdims[4]+1
                across_bdy[1] = -1;
            }
            else
            {
                rdims[4] = ldims[1];
            }
            rdims[1] = rdims[4] - dj;
            if( nj < jextra ) rdims[1]--;
        }
        else
        {
            if( top_j )
            {
                // going across upper periodic bdy in j
                rdims[1]      = gdims[1];
                facedims[4]   = gdims[1];
                across_bdy[1] = 1;
            }
            else
            {
                rdims[1] = ldims[4];
            }
            facedims[1] = facedims[4];
            rdims[4]    = rdims[1] + dj;
            if( nj < jextra ) rdims[4]++;
            if( gperiodic[1] && nj == dijk[1] - 2 ) rdims[4]++;  // +1 because next proc is on periodic bdy
        }
    }
    if( 0 != dijk[2] )
    {
        pto = ( pto + dijk[2] * pijk[1] + np ) % np;
        assert( pto >= 0 && pto < np );
        if( -1 == dijk[2] )
        {
            facedims[5] = facedims[2];
            rdims[5]    = ldims[2];
            rdims[2] -= dk;
            if( pto / pijk[1] < kextra ) rdims[2]--;
        }
        else
        {
            facedims[2] = facedims[5];
            rdims[2]    = ldims[5];
            rdims[5] += dk;
            if( pto / pijk[1] < kextra ) rdims[5]++;
        }
    }

    assert( -1 == pto || ( rdims[0] >= gdims[0] && rdims[3] <= gdims[3] ) );
    assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && bot_j ) ) ) );
    assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) );
    assert( -1 == pto || ( facedims[0] >= rdims[0] && facedims[3] <= rdims[3] ) );
    assert( -1 == pto ||
            ( ( facedims[1] >= rdims[1] || ( gperiodic[1] && rdims[4] == gdims[4] && facedims[1] == gdims[1] ) ) ) );
    assert( -1 == pto || ( facedims[4] <= rdims[4] ) );
    assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) );
    assert( -1 == pto || ( facedims[0] >= ldims[0] && facedims[3] <= ldims[3] ) );
    assert( -1 == pto || ( facedims[1] >= ldims[1] && facedims[4] <= ldims[4] ) );
    assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) );

    return MB_SUCCESS;
}

ErrorCode ScdInterface::get_neighbor_sqijk( int np,
                                            int pfrom,
                                            const int* const gdims,
                                            const int* const gperiodic,
                                            const int* const dijk,
                                            int& pto,
                                            int* rdims,
                                            int* facedims,
                                            int* across_bdy )
{
    if( gperiodic[0] || gperiodic[1] || gperiodic[2] ) return MB_FAILURE;

    pto           = -1;
    across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;
    int pijk[3], lperiodic[3], ldims[6];
    ErrorCode rval = compute_partition_sqijk( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
    if( MB_SUCCESS != rval ) return rval;
    assert( pijk[0] * pijk[1] * pijk[2] == np );
    pto         = -1;
    bool top[3] = { false, false, false }, bot[3] = { false, false, false };
    // nijk: rank in i/j/k direction
    int nijk[3] = { pfrom % pijk[0], ( pfrom % ( pijk[0] * pijk[1] ) ) / pijk[0], pfrom / ( pijk[0] * pijk[1] ) };

    for( int i = 0; i < 3; i++ )
    {
        if( nijk[i] == pijk[i] - 1 ) top[i] = true;
        if( !nijk[i] ) bot[i] = true;
        if( ( !gperiodic[i] && bot[i] && -1 == dijk[i] ) ||  // downward && not periodic
            ( !gperiodic[i] && top[i] && 1 == dijk[i] ) )    // upward && not periodic
            return MB_SUCCESS;
    }

    std::copy( ldims, ldims + 6, facedims );
    std::copy( ldims, ldims + 6, rdims );
    pto = pfrom;
    int delijk[3], extra[3];
    // nijk_to: rank of pto in i/j/k direction
    int nijk_to[3];
    for( int i = 0; i < 3; i++ )
    {
        delijk[i]  = ( gdims[i + 3] == gdims[i] ? 0 : ( gdims[i + 3] - gdims[i] ) / pijk[i] );
        extra[i]   = ( gdims[i + 3] - gdims[i] ) % delijk[i];
        nijk_to[i] = ( nijk[i] + dijk[i] + pijk[i] ) % pijk[i];
    }
    pto = nijk_to[2] * pijk[0] * pijk[1] + nijk_to[1] * pijk[0] + nijk_to[0];
    assert( pto >= 0 && pto < np );
    for( int i = 0; i < 3; i++ )
    {
        if( 0 != dijk[i] )
        {
            if( -1 == dijk[i] )
            {
                facedims[i + 3] = facedims[i];
                if( bot[i] )
                {
                    // going across lower periodic bdy in i
                    rdims[i + 3]  = gdims[i + 3] + 1;  // +1 because ldims[4] on remote proc is gdims[4]+1
                    across_bdy[i] = -1;
                }
                else
                {
                    rdims[i + 3] = ldims[i];
                }
                rdims[i] = rdims[i + 3] - delijk[i];
                if( nijk[i] < extra[i] ) rdims[i]--;
            }
            else
            {
                if( top[i] )
                {
                    // going across upper periodic bdy in i
                    rdims[i]        = gdims[i];
                    facedims[i + 3] = gdims[i];
                    across_bdy[i]   = 1;
                }
                else
                {
                    rdims[i] = ldims[i + 3];
                }
                facedims[i]  = facedims[i + 3];
                rdims[i + 3] = rdims[i] + delijk[i];
                if( nijk[i] < extra[i] ) rdims[i + 3]++;
                if( gperiodic[i] && nijk[i] == dijk[i] - 2 ) rdims[i + 3]++;  // +1 because next proc is on periodic bdy
            }
        }
    }

    assert( -1 != pto );
#ifndef NDEBUG
    for( int i = 0; i < 3; i++ )
    {
        assert( ( rdims[i] >= gdims[i] && ( rdims[i + 3] <= gdims[i + 3] || ( across_bdy[i] && bot[i] ) ) ) );
        assert( ( ( facedims[i] >= rdims[i] ||
                    ( gperiodic[i] && rdims[i + 3] == gdims[i + 3] && facedims[i] == gdims[i] ) ) ) );
        assert( ( facedims[i] >= ldims[i] && facedims[i + 3] <= ldims[i + 3] ) );
    }
#endif

    return MB_SUCCESS;
}

ErrorCode ScdInterface::get_neighbor_alljorkori( int np,
                                                 int pfrom,
                                                 const int* const gdims,
                                                 const int* const gperiodic,
                                                 const int* const dijk,
                                                 int& pto,
                                                 int* rdims,
                                                 int* facedims,
                                                 int* across_bdy )
{
    ErrorCode rval = MB_SUCCESS;
    pto            = -1;
    if( np == 1 ) return MB_SUCCESS;

    int pijk[3], lperiodic[3], ldims[6];
    rval = compute_partition_alljorkori( np, pfrom, gdims, gperiodic, ldims, lperiodic, pijk );
    if( MB_SUCCESS != rval ) return rval;

    int ind       = -1;
    across_bdy[0] = across_bdy[1] = across_bdy[2] = 0;

    for( int i = 0; i < 3; i++ )
    {
        if( pijk[i] > 1 )
        {
            ind = i;
            break;
        }
    }

    assert( -1 < ind );

    if( !dijk[ind] )
        // no neighbor, pto is already -1, return
        return MB_SUCCESS;

    bool is_periodic = ( ( gperiodic[0] && ind == 0 ) || ( gperiodic[1] && ind == 1 ) );
    if( dijk[( ind + 1 ) % 3] || dijk[( ind + 2 ) % 3] ||                   // stepping in either other two directions
        ( !is_periodic && ldims[ind] == gdims[ind] && dijk[ind] == -1 ) ||  // lower side and going lower
        ( !is_periodic && ldims[3 + ind] >= gdims[3 + ind] &&
          dijk[ind] == 1 ) )  // not >= because ldims is only > gdims when periodic;
                              // higher side and going higher
        return MB_SUCCESS;

    std::copy( ldims, ldims + 6, facedims );
    std::copy( ldims, ldims + 6, rdims );

    int dind  = ( gdims[ind + 3] - gdims[ind] ) / np;
    int extra = ( gdims[ind + 3] - gdims[ind] ) % np;
    if( -1 == dijk[ind] && pfrom )
    {
        // actual left neighbor
        pto               = pfrom - 1;  // no need for %np, because pfrom > 0
        facedims[ind + 3] = facedims[ind];
        rdims[ind + 3]    = ldims[ind];
        rdims[ind]        = rdims[ind + 3] - dind - ( pto < extra ? 1 : 0 );
    }
    else if( 1 == dijk[ind] && pfrom < np - 1 )
    {
        // actual right neighbor
        pto            = pfrom + 1;
        facedims[ind]  = facedims[ind + 3];
        rdims[ind]     = ldims[ind + 3];
        rdims[ind + 3] = rdims[ind] + dind + ( pto < extra ? 1 : 0 );
        if( is_periodic && pfrom == np - 2 ) rdims[ind + 3]++;  // neighbor is on periodic bdy
    }
    else if( -1 == dijk[ind] && !pfrom && gperiodic[ind] )
    {
        // downward across periodic bdy
        pto               = np - 1;
        facedims[ind + 3] = facedims[ind] = gdims[ind];  // by convention, facedims is within gdims, so lower value
        rdims[ind + 3] =
            gdims[ind + 3] + 1;  // by convention, local dims one greater than gdims to indicate global lower value
        rdims[ind]      = rdims[ind + 3] - dind - 1;
        across_bdy[ind] = -1;
    }
    else if( 1 == dijk[ind] && pfrom == np - 1 && is_periodic )
    {
        // right across periodic bdy
        pto               = 0;
        facedims[ind + 3] = facedims[ind] = gdims[ind];  // by convention, facedims is within gdims, so lowest value
        rdims[ind]                        = gdims[ind];
        rdims[ind + 3]                    = rdims[ind] + dind + ( pto < extra ? 1 : 0 );
        across_bdy[ind]                   = 1;
    }

    assert( -1 == pto || ( rdims[0] >= gdims[0] && ( rdims[3] <= gdims[3] || ( across_bdy[0] && !pfrom ) ) ) );
    assert( -1 == pto || ( rdims[1] >= gdims[1] && ( rdims[4] <= gdims[4] || ( across_bdy[1] && !pfrom ) ) ) );
    assert( -1 == pto || ( rdims[2] >= gdims[2] && rdims[5] <= gdims[5] ) );
    assert( -1 == pto || ( facedims[0] >= rdims[0] && facedims[3] <= rdims[3] ) );
    assert( -1 == pto || ( facedims[1] >= rdims[1] && facedims[4] <= rdims[4] ) );
    assert( -1 == pto || ( facedims[2] >= rdims[2] && facedims[5] <= rdims[5] ) );
    assert( -1 == pto || ( facedims[0] >= ldims[0] && facedims[3] <= ldims[3] ) );
    assert( -1 == pto || ( facedims[1] >= ldims[1] && facedims[4] <= ldims[4] ) );
    assert( -1 == pto || ( facedims[2] >= ldims[2] && facedims[5] <= ldims[5] ) );

    return rval;
}

//! get shared vertices for alljorkori partition scheme
#ifndef MOAB_HAVE_MPI
ErrorCode ScdInterface::get_shared_vertices( ParallelComm*,<--- The function 'get_shared_vertices' is never used.
                                             ScdBox*,
                                             std::vector< int >&,
                                             std::vector< int >&,
                                             std::vector< int >& )
{
    return MB_FAILURE;
#else
ErrorCode ScdInterface::get_shared_vertices( ParallelComm* pcomm,
                                             ScdBox* box,
                                             std::vector< int >& procs,
                                             std::vector< int >& offsets,
                                             std::vector< int >& shared_indices )
{
    // get index of partitioned dimension
    const int* ldims = box->box_dims();
    ErrorCode rval;
    int ijkrem[6], ijkface[6], across_bdy[3];

    for( int k = -1; k <= 1; k++ )
    {
        for( int j = -1; j <= 1; j++ )
        {
            for( int i = -1; i <= 1; i++ )
            {
                if( !i && !j && !k ) continue;
                int pto;
                int dijk[] = { i, j, k };
                rval = get_neighbor( pcomm->proc_config().proc_size(), pcomm->proc_config().proc_rank(),
                                     box->par_data(), dijk, pto, ijkrem, ijkface, across_bdy );
                if( MB_SUCCESS != rval ) return rval;
                if( -1 != pto )
                {
                    if( procs.empty() || pto != *procs.rbegin() )
                    {
                        procs.push_back( pto );
                        offsets.push_back( shared_indices.size() );
                    }
                    rval = get_indices( ldims, ijkrem, across_bdy, ijkface, shared_indices );
                    if( MB_SUCCESS != rval ) return rval;

                        // check indices against known #verts on local and remote
                        // begin of this block is shared_indices[*offsets.rbegin()], end is
                        // shared_indices.end(), halfway is
                        // (shared_indices.size()-*offsets.rbegin())/2
#ifndef NDEBUG
                    int start_idx = *offsets.rbegin(), end_idx = shared_indices.size(),
                        mid_idx = ( start_idx + end_idx ) / 2;

                    int num_local_verts = ( ldims[3] - ldims[0] + 1 ) * ( ldims[4] - ldims[1] + 1 ) *
                                          ( -1 == ldims[2] && -1 == ldims[5] ? 1 : ( ldims[5] - ldims[2] + 1 ) ),
                        num_remote_verts = ( ijkrem[3] - ijkrem[0] + 1 ) * ( ijkrem[4] - ijkrem[1] + 1 ) *
                                           ( -1 == ijkrem[2] && -1 == ijkrem[5] ? 1 : ( ijkrem[5] - ijkrem[2] + 1 ) );

                    assert(
                        *std::min_element( &shared_indices[start_idx], &shared_indices[mid_idx] ) >= 0 &&
                        *std::max_element( &shared_indices[start_idx], &shared_indices[mid_idx] ) < num_local_verts &&
                        *std::min_element( &shared_indices[mid_idx], &shared_indices[end_idx] ) >= 0 &&
                        *std::max_element( &shared_indices[mid_idx], &shared_indices[end_idx] ) < num_remote_verts );
#endif
                }
            }
        }
    }

    offsets.push_back( shared_indices.size() );

    return MB_SUCCESS;
#endif
}

std::ostream& operator<<( std::ostream& str, const ScdParData& pd )
{
    str << "Partition method = " << ScdParData::PartitionMethodNames[pd.partMethod] << ", gDims = (" << pd.gDims[0]
        << "," << pd.gDims[1] << "," << pd.gDims[2] << ")-(" << pd.gDims[3] << "," << pd.gDims[4] << "," << pd.gDims[5]
        << "), gPeriodic = (" << pd.gPeriodic[0] << "," << pd.gPeriodic[1] << "," << pd.gPeriodic[2] << "), pDims = ("
        << pd.pDims[0] << "," << pd.pDims[1] << "," << pd.pDims[2] << ")" << std::endl;
    return str;
}

}  // namespace moab