ScdInterface.hpp

Go to the documentation of this file.

#ifndef SCD_INTERFACE
#define SCD_INTERFACE

#include "moab/Interface.hpp"
#include "moab/HomXform.hpp"

#include <iostream>
#include <vector>    
#include "assert.h"

namespace moab {

class StructuredElementSeq;
class EntitySequence;
class ScdVertexData;
class EntitySequence;
class ScdBox;
class ParallelComm;


class ScdParData {
public:
  ScdParData() : partMethod(NOPART), pComm(NULL) {
    gDims[0] = gDims[1] = gDims[2] = gDims[3] = gDims[4] = gDims[5] = 0;
    gPeriodic[0] = gPeriodic[1] = gPeriodic[2] = 0;
    pDims[0] = pDims[1] = pDims[2] = 0;
  }

  enum PartitionMethod {ALLJORKORI = 0, ALLJKBAL, SQIJ, SQJK, SQIJK, TRIVIAL, RCBZOLTAN, NOPART};

  static const char *PartitionMethodNames[NOPART + 1];

  int partMethod;
  
  int gDims[6];

  int gPeriodic[3];

  int pDims[3];

  ParallelComm *pComm;
};
  
class ScdInterface 
{
public:
  friend class ScdBox;


  ScdInterface(Interface *impl, bool find_boxes = false);
  
    // Destructor
  ~ScdInterface();

  Interface *impl() const;
  

  ErrorCode construct_box(HomCoord low, HomCoord high, const double * const coords, unsigned int num_coords,
                          ScdBox *& new_box, int * const lperiodic = NULL, 
                          ScdParData * const par_data = NULL,
                          bool assign_global_ids = false, int resolve_shared_ents = -1);


  ErrorCode create_scd_sequence(HomCoord low, HomCoord high, EntityType type,
                                int starting_id, ScdBox *&new_box, 
                                int *is_periodic = NULL);


  ErrorCode find_boxes(std::vector<ScdBox*> &boxes);


  ErrorCode find_boxes(Range &boxes);


  ErrorCode get_boxes(std::vector<ScdBox*> &boxes);


  Tag box_dims_tag(bool create_if_missing = true);


  Tag global_box_dims_tag(bool create_if_missing = true);


  Tag part_method_tag(bool create_if_missing = true);


  Tag box_periodic_tag(bool create_if_missing = true);


  Tag box_set_tag(bool create_if_missing = true);


  ScdBox *get_scd_box(EntityHandle eh);


  static ErrorCode compute_partition(int np, int nr, const ScdParData &par_data,
                                     int *ldims, int *lperiodic = NULL, int *pdims = NULL);
  

  static ErrorCode get_neighbor(int np, int nr, const ScdParData &spd, const int * const dijk,
                                int &pto, int *rdims, int *facedims, int *across_bdy);
  

  ErrorCode tag_shared_vertices(ParallelComm *pcomm, EntityHandle seth);
  

  ErrorCode tag_shared_vertices(ParallelComm *pcomm, ScdBox *box);
  
protected:
  ErrorCode remove_box(ScdBox *box);
  
  ErrorCode add_box(ScdBox *box);
  
private:

  ErrorCode create_box_set(const HomCoord low, const HomCoord high,
                           EntityHandle &scd_set,
                           int *is_periodic = NULL);
  

  inline static ErrorCode compute_partition_alljorkori(int np, int nr,
                                                       const int * const gijk, const int * const gperiodic, 
                                                       int *lijk, int *lperiodic, int *pijk);
  

  inline static ErrorCode compute_partition_alljkbal(int np, int nr, const int * const gijk, const int * const gperiodic, 
                                                     int *lijk, int *lperiodic, int *pijk);


  inline static ErrorCode compute_partition_sqij(int np, int nr, const int * const gijk, const int * const gperiodic, 
                                                 int *lijk, int *lperiodic, int *pijk);
  

  inline static ErrorCode compute_partition_sqjk(int np, int nr, const int * const gijk, const int * const gperiodic, 
                                                 int *lijk, int *lperiodic, int *pijk);


  inline static ErrorCode compute_partition_sqijk(int np, int nr, const int * const gijk, const int * const gperiodic, 
                                                  int *lijk, int *lperiodic, int *pijk);


  static ErrorCode get_shared_vertices(ParallelComm *pcomm, ScdBox *box, std::vector<int> &procs,
                                       std::vector<int> &offsets, std::vector<int> &shared_indices);

  static ErrorCode get_indices(const int * const ldims, const int * const rdims, const int * const across_bdy, 
                               int *face_dims, std::vector<int> &shared_indices);
  
  static ErrorCode 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);
  
  static ErrorCode 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);
  
  static ErrorCode 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);
  
  static ErrorCode 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);
  
  static ErrorCode 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);
  
  static int gtol(const int *gijk, int i, int j, int k);

  ErrorCode assign_global_ids(ScdBox *box);
  
  Interface *mbImpl;

  bool searchedBoxes;
  
  std::vector<ScdBox*> scdBoxes;

  Tag boxPeriodicTag;

  Tag boxDimsTag;

  Tag globalBoxDimsTag;

  Tag partMethodTag;

  Tag boxSetTag;
  
};

class ScdBox 
{
  friend class ScdInterface;
  
public:

  ~ScdBox();

  inline ScdInterface *sc_impl() const;

    /* Add a vertex box to the element sequence referenced by this box.  The passed in vbox must
     * be a vertex box, with parametric extents no larger than that of this box.  This vbox is
     * oriented to this box by matching parameters from1-from3 in vbox to to1-to3 in this box.
     * If bb_input is true, only the part of the vertex sequence between bb_min and bb_max is referenced
     * \param vbox The vertex box being added to this box
     * \param from1 1st reference point on vbox
     * \param to1 1st reference point on this box
     * \param from2 2nd reference point on vbox
     * \param to2 2nd reference point on this box
     * \param from3 3rd reference point on vbox
     * \param to3 3rd reference point on this box
     * \param bb_input If true, subsequent parameters list extents of vbox referenced
     * \param bb_min Lower corner of rectangle referenced
     * \param bb_max Upper corner of rectangle referenced
     */
  ErrorCode add_vbox(ScdBox *vbox,
                     HomCoord from1, HomCoord to1, 
                     HomCoord from2, HomCoord to2,
                     HomCoord from3, HomCoord to3,
                     bool bb_input = false,
                     const HomCoord &bb_min = HomCoord::unitv[0],
                     const HomCoord &bb_max = HomCoord::unitv[0]);


  bool boundary_complete() const;

  inline int box_dimension() const;
  
  inline EntityHandle start_vertex() const;
  

  inline EntityHandle start_element() const;

    /* Number of elements is (boxSize[0]-1)(boxSize[1]-1)(boxSize[2]-1)
     */
  inline int num_elements() const;
  
    /* Number of vertices is boxSize[0] * boxSize[1] * boxSize[2]
     */
  inline int num_vertices() const;
  

  inline const int *box_dims() const;
  
  inline HomCoord box_min() const;
  
  inline HomCoord box_max() const;
  
  inline HomCoord box_size() const;
  

  inline void box_size(int *ijk) const;
  

  void box_size(int &i, int &j, int &k) const;
  

  EntityHandle get_element(HomCoord ijk) const;
  

  EntityHandle get_element(int i, int j = 0, int k = 0) const;
  

  EntityHandle get_vertex(HomCoord ijk) const;
  

  EntityHandle get_vertex(int i, int j = 0, int k = 0) const;
  

  ErrorCode get_params(EntityHandle ent, int &i, int &j, int &k, int &dir) const;
  

  ErrorCode get_params(EntityHandle ent, int &i, int &j, int &k) const;
  

  ErrorCode get_params(EntityHandle ent, HomCoord &ijkd) const;
  
  ErrorCode get_adj_edge_or_face(int dim, int i, int j, int k, int dir, EntityHandle &ent,
                                 bool create_if_missing = true) const;


  bool contains(int i, int j, int k) const;


  bool contains(const HomCoord ijk) const;

  void box_set(EntityHandle this_set);
  EntityHandle box_set();


  ErrorCode get_coordinate_arrays(double *&xc, double *&yc, double *&zc);
  

  ErrorCode get_coordinate_arrays(const double *&xc, const double *&yc, const double *&zc) const;


  bool locally_periodic_i() const;
  

  bool locally_periodic_j() const;
  

  bool locally_periodic_k() const;
  

  void locally_periodic(bool lperiodic[3]);


  const int *locally_periodic() const;
 

  ScdParData &par_data() {return parData;}
  

  const ScdParData &par_data() const {return parData;}
  

  void par_data(const ScdParData &par_datap) {parData = par_datap;}
  
private:

  ScdBox(ScdInterface *sc_impl, EntityHandle box_set,
         EntitySequence *seq1, EntitySequence *seq2 = NULL);


  EntityHandle get_vertex_from_seq(int i, int j, int k) const;

  ErrorCode vert_dat(ScdVertexData *vert_dat);
  
  ScdVertexData *vert_dat() const;
  
  ErrorCode elem_seq(EntitySequence *elem_seq);
  
  StructuredElementSeq *elem_seq() const;
  
  void start_vertex(EntityHandle startv);
  
  void start_element(EntityHandle starte);

  ScdInterface *scImpl;

  EntityHandle boxSet;

  ScdVertexData *vertDat;

  StructuredElementSeq *elemSeq;
  
  EntityHandle startVertex;
  
  EntityHandle startElem;

  int boxDims[6];

  int locallyPeriodic[3];

  ScdParData parData;
  
  HomCoord boxSize;

  int boxSizeIJ;
  int boxSizeIJM1;
  int boxSizeIM1;
  
};

inline ErrorCode ScdInterface::compute_partition(int np, int nr, const ScdParData &par_data,
                                     int *ldims, int *lperiodic, int *pdims)
{
  ErrorCode rval = MB_SUCCESS;
  switch (par_data.partMethod) {
    case ScdParData::ALLJORKORI:
    case -1:
        rval = compute_partition_alljorkori(np, nr, par_data.gDims, par_data.gPeriodic, 
                                            ldims, lperiodic, pdims);
        break;
    case ScdParData::ALLJKBAL:
        rval = compute_partition_alljkbal(np, nr, par_data.gDims, par_data.gPeriodic, 
                                          ldims, lperiodic, pdims);
        break;
    case ScdParData::SQIJ:
        rval = compute_partition_sqij(np, nr, par_data.gDims, par_data.gPeriodic, 
                                      ldims, lperiodic, pdims);
        break;
    case ScdParData::SQJK:
        rval = compute_partition_sqjk(np, nr, par_data.gDims, par_data.gPeriodic, 
                                      ldims, lperiodic, pdims);
        break;
    case ScdParData::SQIJK:
        rval = compute_partition_sqijk(np, nr, par_data.gDims, par_data.gPeriodic, 
                                      ldims, lperiodic, pdims);
        break;
    default:
        rval = MB_FAILURE;
        break;
  }

  return rval;
}

inline ErrorCode ScdInterface::compute_partition_alljorkori(int np, int nr,
                                                            const int * const gijk, const int * const gperiodic, 
                                                            int *ldims, int *lperiodic, int *pijk)
{
    // partition *the elements* over the parametric space; 1d partition for now, in the j, k, or i
    // parameters
  int tmp_lp[3], tmp_pijk[3];
  if (!lperiodic) lperiodic = tmp_lp;
  if (!pijk) pijk = tmp_pijk;
  
  for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];
  
  if (np == 1) {
    if (ldims) {
      ldims[0] = gijk[0];
      ldims[3] = gijk[3];
      ldims[1] = gijk[1];
      ldims[4] = gijk[4];
      ldims[2] = gijk[2];
      ldims[5] = gijk[5];
    }
    pijk[0] = pijk[1] = pijk[2] = 1;
  }
  else {
    if (gijk[4] - gijk[1] > np) {
        // partition j over procs
      int dj = (gijk[4] - gijk[1]) / np;
      int extra = (gijk[4] - gijk[1]) % np;
      ldims[1] = gijk[1] + nr*dj + 
          std::min(nr, extra);
      ldims[4] = ldims[1] + dj + (nr < extra ? 1 : 0);

      if (gperiodic[1] && np > 1) {
        lperiodic[1] = 0;
        ldims[4]++;
      }
      
      ldims[2] = gijk[2]; ldims[5] = gijk[5];
      ldims[0] = gijk[0]; ldims[3] = gijk[3];
      pijk[0] = pijk[2] = 1;
      pijk[1] = np;
    }
    else if (gijk[5] - gijk[2] > np) {
        // partition k over procs
      int dk = (gijk[5] - gijk[2]) / np;
      int extra = (gijk[5] - gijk[2]) % np;
      ldims[2] = gijk[2] + nr*dk + 
          std::min(nr, extra);
      ldims[5] = ldims[2] + dk + (nr < extra ? 1 : 0);

      ldims[1] = gijk[1]; ldims[4] = gijk[4];
      ldims[0] = gijk[0]; ldims[3] = gijk[3];
      pijk[0] = pijk[1] = 1;
      pijk[2] = np;
    }
    else if (gijk[3] - gijk[0] > np) {
        // partition i over procs
      int di = (gijk[3] - gijk[0]) / np;
      int extra = (gijk[3] - gijk[0]) % np;
      ldims[0] = gijk[0] + nr*di + 
          std::min(nr, extra);
      ldims[3] = ldims[0] + di + (nr < extra ? 1 : 0);

      if (gperiodic[0] && np > 1) {
        lperiodic[0] = 0;
        ldims[3]++;
      }

      ldims[2] = gijk[2]; ldims[5] = gijk[5];
      ldims[1] = gijk[1]; ldims[4] = gijk[4];

      pijk[1] = pijk[2] = 1;
      pijk[0] = np;
    }
    else {
        // Couldn't find a suitable partition...
      return MB_FAILURE;
    }
  }
  
  return MB_SUCCESS;
}

inline ErrorCode ScdInterface::compute_partition_alljkbal(int np, int nr,
                                                          const int * const gijk, const int * const gperiodic, 
                                                          int *ldims, int *lperiodic, int *pijk)
{
  int tmp_lp[3], tmp_pijk[3];
  if (!lperiodic) lperiodic = tmp_lp;
  if (!pijk) pijk = tmp_pijk;

  for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];

  if (np == 1) {
    if (ldims) {
      ldims[0] = gijk[0];
      ldims[3] = gijk[3];
      ldims[1] = gijk[1];
      ldims[4] = gijk[4];
      ldims[2] = gijk[2];
      ldims[5] = gijk[5];
    }
    pijk[0] = pijk[1] = pijk[2] = 1;
  }
  else {
      // improved, possibly 2-d partition
    std::vector<double> kfactors;
    kfactors.push_back(1);
    int K = gijk[5] - gijk[2];
    for (int i = 2; i < K; i++) 
      if (!(K%i) && !(np%i)) kfactors.push_back(i);
    kfactors.push_back(K);
  
      // compute the ideal nj and nk
    int J = gijk[4] - gijk[1];
    double njideal = sqrt(((double)(np*J))/((double)K));
    double nkideal = (njideal*K)/J;
  
    int nk, nj;
    if (nkideal < 1.0) {
      nk = 1;
      nj = np;
    }
    else {
      std::vector<double>::iterator vit = std::lower_bound(kfactors.begin(), kfactors.end(), nkideal);
      if (vit == kfactors.begin()) nk = 1;
      else nk = (int)*(--vit);
      nj = np / nk;
    }

    int dk = K / nk;
    int dj = J / nj;
  
    ldims[2] = gijk[2] + (nr % nk) * dk;
    ldims[5] = ldims[2] + dk;
  
    int extra = J % nj;
  
    ldims[1] = gijk[1] + (nr / nk) * dj + std::min(nr / nk, extra);
    ldims[4] = ldims[1] + dj + (nr / nk < extra ? 1 : 0);

    ldims[0] = gijk[0];
    ldims[3] = gijk[3];

    if (gperiodic[1] && np > 1) {
      lperiodic[1] = 0;
      if (nr/nk == nj-1) {
        ldims[1]++;
      }
    }

    pijk[0] = 1; pijk[1] = nj; pijk[2] = nk;
  }
  
  return MB_SUCCESS;
}

inline ErrorCode ScdInterface::compute_partition_sqij(int np, int nr,
                                                      const int * const gijk, const int * const gperiodic, 
                                                      int *ldims, int *lperiodic, int *pijk)
{
  int tmp_lp[3], tmp_pijk[3];
  if (!lperiodic) lperiodic = tmp_lp;
  if (!pijk) pijk = tmp_pijk;

    // square IxJ partition

  for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];
  
  if (np == 1) {
    if (ldims) {
      ldims[0] = gijk[0];
      ldims[3] = gijk[3];
      ldims[1] = gijk[1];
      ldims[4] = gijk[4];
      ldims[2] = gijk[2];
      ldims[5] = gijk[5];
    }
    pijk[0] = pijk[1] = pijk[2] = 1;
  }
  else {
    std::vector<double> pfactors, ppfactors;
    for (int i = 2; i <= np/2; i++) 
      if (!(np%i)) {
        pfactors.push_back(i);
        ppfactors.push_back(((double)(i*i))/np);
      }
    pfactors.push_back(np);
    ppfactors.push_back( (double)np);
    
      // ideally, Px/Py = I/J
    double ijratio = ((double)(gijk[3]-gijk[0]))/((double)(gijk[4]-gijk[1]));
    
    unsigned int ind = std::lower_bound(ppfactors.begin(), ppfactors.end(), ijratio) - ppfactors.begin();
    if (ind && fabs(ppfactors[ind-1]-ijratio) < fabs(ppfactors[ind]-ijratio)) ind--;
    

      // VARIABLES DESCRIBING THE MESH:
      // pi, pj = # procs in i and j directions
      // nri, nrj = my proc's position in i, j directions
      // I, J = # edges/elements in i, j directions
      // iextra, jextra = # procs having extra edge in i/j direction
      // top_i, top_j = if true, I'm the last proc in the i/j direction
      // i, j = # edges locally in i/j direction, *not* including one for iextra/jextra
    int pi = pfactors[ind];
    int pj = np / pi;
    
    int I = (gijk[3] - gijk[0]), J = (gijk[4] - gijk[1]);
    int iextra = I%pi, jextra = J%pj, i = I/pi, j = J/pj;
    int nri = nr % pi, nrj = nr / pi;

    if (ldims) {
      ldims[0] = gijk[0] + i*nri + std::min(iextra, nri);
      ldims[3] = ldims[0] + i + (nri < iextra ? 1 : 0);
      ldims[1] = gijk[1] + j*nrj + std::min(jextra, nrj);
      ldims[4] = ldims[1] + j + (nrj < jextra ? 1 : 0);
      
      ldims[2] = gijk[2];
      ldims[5] = gijk[5];

      if (gperiodic[0] && pi > 1) {
        lperiodic[0] = 0;
        if (nri == pi-1) ldims[3]++;
      }
      if (gperiodic[1] && pj > 1) {
        lperiodic[1] = 0;
        if (nrj == pj-1) ldims[4]++;
      }
      
    }

    pijk[0] = pi; pijk[1] = pj; pijk[2] = 1;
  }  

  return MB_SUCCESS;
}

inline ErrorCode ScdInterface::compute_partition_sqjk(int np, int nr,
                                                      const int * const gijk, const int * const gperiodic, 
                                                      int *ldims, int *lperiodic, int *pijk)
{
  int tmp_lp[3], tmp_pijk[3];
  if (!lperiodic) lperiodic = tmp_lp;
  if (!pijk) pijk = tmp_pijk;

    // square JxK partition
  for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];

  if (np == 1) {
    if (ldims) {
      ldims[0] = gijk[0];
      ldims[3] = gijk[3];
      ldims[1] = gijk[1];
      ldims[4] = gijk[4];
      ldims[2] = gijk[2];
      ldims[5] = gijk[5];
    }
    pijk[0] = pijk[1] = pijk[2] = 1;
  }
  else {
    std::vector<double> pfactors, ppfactors;
    for (int p = 2; p <= np; p++) 
      if (!(np%p)) {
        pfactors.push_back(p);
        ppfactors.push_back(((double)(p*p))/np);
      }
  
      // ideally, Pj/Pk = J/K
    int pj, pk;
    if (gijk[5] == gijk[2]) {
      pk = 1;
      pj = np;
    }
    else {
      double jkratio = ((double)(gijk[4]-gijk[1]))/((double)(gijk[5]-gijk[2]));

      std::vector<double>::iterator vit  = std::lower_bound(ppfactors.begin(), ppfactors.end(), jkratio);
      if (vit == ppfactors.end()) vit--;
      else if (vit != ppfactors.begin() && fabs(*(vit-1)-jkratio) < fabs((*vit)-jkratio)) vit--;
      int ind = vit - ppfactors.begin();
  
      pj = 1;
      if (ind >= 0 && !pfactors.empty()) pfactors[ind];
      pk = np / pj;
    }

    int K = (gijk[5] - gijk[2]), J = (gijk[4] - gijk[1]);
    int jextra = J%pj, kextra = K%pk, j = J/pj, k = K/pk;
    int nrj = nr % pj, nrk = nr / pj;
    ldims[1] = gijk[1] + j*nrj + std::min(jextra, nrj);
    ldims[4] = ldims[1] + j + (nrj < jextra ? 1 : 0);
    ldims[2] = gijk[2] + k*nrk + std::min(kextra, nrk);
    ldims[5] = ldims[2] + k + (nrk < kextra ? 1 : 0);

    ldims[0] = gijk[0];
    ldims[3] = gijk[3];

    if (gperiodic[1] && pj > 1) {
      lperiodic[1] = 0;
      if (nrj == pj-1) ldims[4]++;
    }

    pijk[0] = 1; pijk[1] = pj; pijk[2] = pk;
  }
  
  return MB_SUCCESS;
}

inline ErrorCode ScdInterface::compute_partition_sqijk(int np, int nr,
                                                       const int * const gijk, const int * const gperiodic, 
                                                       int *ldims, int *lperiodic, int *pijk)
{
  if (gperiodic[0] || gperiodic[1] || gperiodic[2]) return MB_FAILURE;
  
  int tmp_lp[3], tmp_pijk[3];
  if (!lperiodic) lperiodic = tmp_lp;
  if (!pijk) pijk = tmp_pijk;

    // square IxJxK partition

  for (int i = 0; i < 3; i++) lperiodic[i] = gperiodic[i];
  
  if (np == 1) {
    if (ldims) 
      for (int i = 0; i < 6; i++) ldims[i] = gijk[i];
    pijk[0] = pijk[1] = pijk[2] = 1;
    return MB_SUCCESS;
  }

  std::vector<int> pfactors;
  pfactors.push_back(1);
  for (int i = 2; i <= np/2; i++) 
    if (!(np%i)) pfactors.push_back(i);
  pfactors.push_back(np);

    // test for IJ, JK, IK
  int IJK[3], dIJK[3];
  for (int i = 0; i < 3; i++)
    IJK[i] = std::max(gijk[3+i] - gijk[i], 1);
    // order IJK from lo to hi
  int lo = 0, hi = 0;
  for (int i = 1; i < 3; i++) {
    if (IJK[i] < IJK[lo]) lo = i;
    if (IJK[i] > IJK[hi]) hi = i;
  }
  if (lo == hi) hi = (lo+1)%3;
  int mid = 3 - lo - hi;
    // search for perfect subdivision of np that balances #cells
  int perfa_best = -1, perfb_best = -1;
  double ratio = 0.0;
  for (int po = 0; po < (int)pfactors.size(); po++) {
    for (int pi = po; pi < (int)pfactors.size() && np/(pfactors[po]*pfactors[pi]) >= pfactors[pi]; pi++) {
      int p3 = std::find(pfactors.begin(), pfactors.end(), np/(pfactors[po]*pfactors[pi])) - pfactors.begin();
      if (p3 == (int)pfactors.size() || pfactors[po]*pfactors[pi]*pfactors[p3] != np) continue; // po*pi should exactly factor np
      assert(po <= pi && pi <= p3);
        // by definition, po <= pi <= p3
      double minl = std::min(std::min(IJK[lo]/pfactors[po], IJK[mid]/pfactors[pi]), IJK[hi]/pfactors[p3]),
          maxl = std::max(std::max(IJK[lo]/pfactors[po], IJK[mid]/pfactors[pi]), IJK[hi]/pfactors[p3]);
      if (minl/maxl > ratio) {
        ratio = minl/maxl;
        perfa_best = po; perfb_best = pi;
      }
    }
  }
  if (perfa_best == -1 || perfb_best == -1) 
    return MB_FAILURE;

    // VARIABLES DESCRIBING THE MESH:
    // pijk[i] = # procs in direction i
    // numr[i] = my proc's position in direction i
    // dIJK[i] = # edges/elements in direction i
    // extra[i]= # procs having extra edge in direction i
    // top[i] = if true, I'm the last proc in direction i
    
  pijk[lo] = pfactors[perfa_best]; pijk[mid] = pfactors[perfb_best]; 
  pijk[hi] = (np/(pfactors[perfa_best]*pfactors[perfb_best]));
  int extra[3] = {0, 0, 0}, numr[3];
  for (int i = 0; i < 3; i++) {
    dIJK[i] = IJK[i] / pijk[i];
    extra[i] = IJK[i]%pijk[i];
  }
  numr[2] = nr / (pijk[0]*pijk[1]);
  int rem = nr % (pijk[0] * pijk[1]);
  numr[1] = rem / pijk[0];
  numr[0] = rem % pijk[0];
  for (int i = 0; i < 3; i++) {
    extra[i] = IJK[i] % dIJK[i];
    ldims[i] = gijk[i] + numr[i]*dIJK[i] + std::min(extra[i], numr[i]);
    ldims[3+i] = ldims[i] + dIJK[i] + (numr[i] < extra[i] ? 1 : 0);
  }

  return MB_SUCCESS;
}

inline int ScdInterface::gtol(const int *gijk, int i, int j, int k) 
{
  return ((k-gijk[2])*(gijk[3]-gijk[0]+1)*(gijk[4]-gijk[1]+1) + (j-gijk[1])*(gijk[3]-gijk[0]+1) + i-gijk[0]);
}

inline ErrorCode ScdInterface::get_indices(const int * const ldims, const int * const rdims, const int * const across_bdy, 
                                           int *face_dims, std::vector<int> &shared_indices) 
{
    // check for going across periodic bdy and face_dims not in my ldims (I'll always be on top in that case)...
  if (across_bdy[0] > 0 && face_dims[0] != ldims[3]) face_dims[0] = face_dims[3] = ldims[3];
  else if (across_bdy[0] < 0 && face_dims[0] != ldims[0]) face_dims[0] = face_dims[3] = ldims[0];
  if (across_bdy[1] > 0 && face_dims[1] != ldims[4]) face_dims[1] = face_dims[4] = ldims[4];
  else if (across_bdy[1] < 0 && face_dims[1] != ldims[1]) face_dims[0] = face_dims[3] = ldims[1];
  
  for (int k = face_dims[2]; k <= face_dims[5]; k++)
    for (int j = face_dims[1]; j <= face_dims[4]; j++)
      for (int i = face_dims[0]; i <= face_dims[3]; i++)
        shared_indices.push_back(gtol(ldims, i, j, k));

  if (across_bdy[0] > 0 && face_dims[0] != rdims[0]) face_dims[0] = face_dims[3] = rdims[0];
  else if (across_bdy[0] < 0 && face_dims[0] != rdims[3]) face_dims[0] = face_dims[3] = rdims[3];
  if (across_bdy[1] > 0 && face_dims[1] != rdims[1]) face_dims[1] = face_dims[4] = rdims[1];
  else if (across_bdy[1] < 0 && face_dims[1] != rdims[4]) face_dims[0] = face_dims[3] = rdims[4];
  
  for (int k = face_dims[2]; k <= face_dims[5]; k++)
    for (int j = face_dims[1]; j <= face_dims[4]; j++)
      for (int i = face_dims[0]; i <= face_dims[3]; i++)
        shared_indices.push_back(gtol(rdims, i, j, k));

  return MB_SUCCESS;
}
  
inline ErrorCode ScdInterface::get_neighbor(int np, int pfrom, const ScdParData &spd, const int * const dijk, 
                                            int &pto, int *rdims, int *facedims, int *across_bdy)
{
  if (!dijk[0] && !dijk[1] && !dijk[2]) {
    // not going anywhere, return
    pto = -1;
    return MB_SUCCESS;
  }
  
  switch (spd.partMethod) {
    case ScdParData::ALLJORKORI:
    case -1:
        return get_neighbor_alljorkori(np, pfrom, spd.gDims, spd.gPeriodic, dijk,
                                       pto, rdims, facedims, across_bdy);
    case ScdParData::ALLJKBAL:
        return get_neighbor_alljkbal(np, pfrom, spd.gDims, spd.gPeriodic, dijk,
                                     pto, rdims, facedims, across_bdy);
    case ScdParData::SQIJ:
        return get_neighbor_sqij(np, pfrom, spd.gDims, spd.gPeriodic, dijk,
                                 pto, rdims, facedims, across_bdy);
    case ScdParData::SQJK:
        return get_neighbor_sqjk(np, pfrom, spd.gDims, spd.gPeriodic, dijk,
                                 pto, rdims, facedims, across_bdy);
    case ScdParData::SQIJK:
        return get_neighbor_sqijk(np, pfrom, spd.gDims, spd.gPeriodic, dijk,
                                  pto, rdims, facedims, across_bdy);
    default:
        break;
  }

  return MB_FAILURE;
}

inline ErrorCode ScdInterface::tag_shared_vertices(ParallelComm *pcomm, EntityHandle seth) 
{
  ScdBox *box = get_scd_box(seth);
  if (!box) {
      // look for contained boxes
    Range tmp_range;
    ErrorCode rval = mbImpl->get_entities_by_type(seth, MBENTITYSET, tmp_range);
    if (MB_SUCCESS != rval) return rval;
    for (Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); rit++) {
      box = get_scd_box(*rit);
      if (box) break;
    }
  }
  
  if (!box) return MB_FAILURE;

  return tag_shared_vertices(pcomm, box);
}

inline ScdInterface *ScdBox::sc_impl() const 
{
  return scImpl;
}

inline EntityHandle ScdBox::start_vertex() const
{
  return startVertex;
}
    
inline void ScdBox::start_vertex(EntityHandle startv) 
{
  startVertex = startv;
}
    
inline EntityHandle ScdBox::start_element() const
{
  return startElem;
}
    
inline void ScdBox::start_element(EntityHandle starte) 
{
  startElem = starte;
}
    
inline int ScdBox::num_elements() const
{
  return (!startElem ? 0 : 
          (boxSize[0]- (locallyPeriodic[0] ? 0 : 1)) * 
          (-1 == boxSize[1] ? 1 : (boxSize[1]-(locallyPeriodic[1] ? 0 : 1))) * 
          (boxSize[2] == -1 || boxSize[2] == 1 ? 1 : (boxSize[2]-(locallyPeriodic[2] ? 0 : 1))));
}
    
inline int ScdBox::num_vertices() const
{
  return boxSize[0] * (!boxSize[1] ? 1 : boxSize[1]) * 
      (!boxSize[2] ? 1 : boxSize[2]);
}
    
inline const int *ScdBox::box_dims() const 
{
  return boxDims;
}

inline HomCoord ScdBox::box_min() const 
{
  return HomCoord(boxDims, 3);
}

inline HomCoord ScdBox::box_max() const 
{
  return HomCoord(boxDims+3, 3);
}

inline HomCoord ScdBox::box_size() const 
{
  return boxSize;
}

inline void ScdBox::box_size(int *ijk) const 
{
  ijk[0] = boxSize[0];
  ijk[1] = boxSize[1];
  ijk[2] = boxSize[2];
}

inline void ScdBox::box_size(int &i, int &j, int &k) const 
{
  i = boxSize[0];
  j = boxSize[1];
  k = boxSize[2];
}
  
inline EntityHandle ScdBox::get_element(int i, int j, int k) const 
{
  return (!startElem ? 0 : 
          startElem + (k-boxDims[2])*boxSizeIJM1 + (j-boxDims[1])*boxSizeIM1 + i-boxDims[0]);
}

inline EntityHandle ScdBox::get_element(HomCoord ijk) const 
{
  return get_element(ijk[0], ijk[1], ijk[2]);
}
  
inline EntityHandle ScdBox::get_vertex(int i, int j, int k) const 
{
  return (vertDat ? startVertex + (boxDims[2] == -1 && boxDims[5] == -1 ? 0 : (k-boxDims[2]))*boxSizeIJ +
      (boxDims[1] == -1 && boxDims[4] == -1 ? 0 : (j-boxDims[1]))*boxSize[0] + i-boxDims[0] : get_vertex_from_seq(i, j, k));
}

inline EntityHandle ScdBox::get_vertex(HomCoord ijk) const 
{
  return get_vertex(ijk[0], ijk[1], ijk[2]);
}
  
inline bool ScdBox::contains(const HomCoord ijk) const
{
  return (ijk >= HomCoord(boxDims, 3) && 
          ijk <= HomCoord(boxDims+3, 3));
}

inline bool ScdBox::contains(int i, int j, int k) const
{
  return contains(HomCoord(i, j, k));
}

inline void ScdBox::box_set(EntityHandle this_set) 
{
  boxSet = this_set;
}

inline EntityHandle ScdBox::box_set() 
{
  return boxSet;
}
    
inline ScdVertexData *ScdBox::vert_dat() const
{
  return vertDat;
}
  
inline StructuredElementSeq *ScdBox::elem_seq() const
{
  return elemSeq;
}
  
inline ErrorCode ScdBox::get_params(EntityHandle ent, int &i, int &j, int &k, int &dir) const
{
  HomCoord hc;
  ErrorCode rval = get_params(ent, hc);
  if (MB_SUCCESS == rval) {
    i = hc[0];
    j = hc[1];
    k = hc[2];
    dir = hc[3];
  }
  
  return rval;
}

inline ErrorCode ScdBox::get_params(EntityHandle ent, int &i, int &j, int &k) const 
{
  HomCoord hc;
  ErrorCode rval = get_params(ent, hc);
  if (MB_SUCCESS == rval) {
    i = hc[0];
    j = hc[1];
    k = hc[2];
  }
  
  return rval;
}

inline bool ScdBox::locally_periodic_i() const 
{
  return locallyPeriodic[0];
}

inline bool ScdBox::locally_periodic_j() const 
{
  return locallyPeriodic[1];
}

inline bool ScdBox::locally_periodic_k() const 
{
  return locallyPeriodic[2];
}

inline void ScdBox::locally_periodic(bool lperiodic[3])
{
   for (int i = 0; i < 3; i++) 
    locallyPeriodic[i] = lperiodic[i];
}

inline const int *ScdBox::locally_periodic() const
{
  return locallyPeriodic;
}
 
inline 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
#endif