ScdInterface.cpp

Go to the documentation of this file.

#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 = nullptr;
            // std::remove_if(scdBoxes.begin(), scdBoxes.end(),
            // std::bind2nd(std::equal_to<ScdBox*>(), dum) ) ;
            std::remove_if( scdBoxes.begin(), scdBoxes.end(),
                            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 )
{
    // 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 )
{
    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 ) &&
                    ( elemSeq->sdata()->max_params() + HomCoord( 1, 1, 1 ) ) == HomCoord( boxDims, 3 ) );
        }
        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 )
{
    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*,
                                             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