parallel element-to-element orientation arrays

orientation.cc

#include <apfMesh.h>
#include <apfNumbering.h>
#include <PCU.h>

static int getOrientationCode(
    int side_i,
    long* verts,
    int other_side_i,
    long* other_verts)
{
  return 42;
}

static void sendOrientation(
    apf::Mesh* m,
    apf::MeshEntity* side,
    int side_i,
    long* gids,
    int nverts)
{
  apf::Copy other = apf::getOtherCopy(m, side);
  PCU_COMM_PACK(other.peer, other.entity);
  PCU_COMM_PACK(other.peer, side_i);
  for (int i = 0; i < nverts; ++i)
    PCU_COMM_PACK(other.peer, gids[i]);
}

static void receiveOrientation(
    apf::GlobalNumbering* gn,
    apf::Numbering* n,
    int* o)
{
  apf::Mesh* m = apf::getMesh(n);
  apf::MeshEntity* side;
  PCU_COMM_UNPACK(side);
  int other_side_i;
  PCU_COMM_UNPACK(other_side_i);
  apf::MeshEntity* element = m->getUpward(side, 0);
  apf::Downward sides;
  int dim = m->getDimension();
  int nsides = m->getDownward(element, dim - 1, sides);
  int side_i = apf::findIn(sides, nsides, side);
  apf::NewArray<long> gids;
  apf::getElementNumbers(gn, element, gids);
  int nverts = apf::Mesh::adjacentCount[m->getType(element)][0];
  apf::NewArray<long> other_gids(nverts);
  for (int i = 0; i < nverts; ++i)
    PCU_COMM_UNPACK(other_gids[i]);
  int lid = apf::getNumber(n, element, 0, 0);
  o[lid * nsides + side_i] = getOrientationCode(
      side_i, &gids[0], other_side_i, &other_gids[0]);
}

static int getLocalOrientation(
    apf::GlobalNumbering* gn,
    apf::Numbering* n,
    apf::MeshEntity* element,
    apf::MeshEntity* side,
    int side_i,
    long* gids)
{
  apf::Mesh* m = apf::getMesh(gn);
  apf::Up elements;
  m->getUp(side, elements);
  if (elements.n == 1)
    return -1;
  int i = apf::findIn(elements.e, 2, element);
  apf::MeshEntity* other_element = elements.e[1 - i];
  apf::Downward other_sides;
  int dim = m->getDimension();
  int nsides = m->getDownward(other_element, dim - 1, other_sides);
  int other_side_i = apf::findIn(other_sides, nsides, side);
  apf::NewArray<long> other_gids;
  apf::getElementNumbers(gn, other_element, other_gids);
  return getOrientationCode(
      side_i, gids, other_side_i, &other_gids[0]);
}

static void getOrientations(
    apf::GlobalNumbering* gn,
    apf::Numbering* n,
    int* o,
    apf::MeshEntity* e)
{
  apf::Mesh* m = apf::getMesh(n);
  apf::Downward sides;
  int dim = m->getDimension();
  int nsides = m->getDownward(e, dim - 1, sides);
  apf::NewArray<long> gids;
  apf::getElementNumbers(gn, e, gids);
  int lid = apf::getNumber(n, e, 0, 0);
  int nverts = apf::Mesh::adjacentCount[m->getType(e)][0];
  for (int i = 0; i < nsides; ++i) {
    apf::MeshEntity* side = sides[i];
    if (m->isShared(side))
      sendOrientation(m, side, i, &gids[0], nverts);
    else
      o[lid * nsides + i] = getLocalOrientation(
          gn, n, e, side, i, &gids[0]);
  }
}

int* getInterElementOrientations(apf::Mesh* m)
{
  int dim = m->getDimension();
  int type = apf::getFirstType(m, dim);
  int nsides = apf::Mesh::adjacentCount[type][dim - 1];
  int* o = new int[m->count(dim) * nsides];
  apf::GlobalNumbering* gn = apf::makeGlobal(
      apf::numberOwnedNodes(m, "ieo_node"));
  apf::synchronize(gn);
  apf::Numbering* n = apf::numberElements(m, "ieo_elem");
  PCU_Comm_Begin();
  apf::MeshIterator* it = m->begin(dim);
  apf::MeshEntity* e;
  while ((e = m->iterate(it)))
    getOrientations(gn, n, o, e);
  m->end(it);
  PCU_Comm_Send();
  while (PCU_Comm_Receive())
    receiveOrientation(gn, n, o);
  apf::destroyNumbering(n);
  apf::destroyGlobalNumbering(gn);
  return o;
}