#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: tri2dMovement.c,v 1.4 2000/01/10 03:54:23 knepley Exp $";
#endif

/* Mesh movement for 2d triangular grids */
#include "src/gvec/meshMovementImpl.h"           /*I "meshMovement.h" I*/
#include "src/mesh/impls/triangular/2d/2dimpl.h" /*I "mesh.h" I*/
#include "src/gvec/gvecimpl.h"                   /* For MeshMoveMesh() */
#include "tri2dMovement.h"

#undef  __FUNCT__
#define __FUNCT__ "MeshMoverSetup_Triangular_2D"
int MeshMoverSetup_Triangular_2D(MeshMover mover) {
  Mesh               mesh       = mover->mesh;
  Mesh               movingMesh = mesh;
  int                dim        = mesh->dim;
  int                numCorners = mesh->numCorners;
  DiscretizationType dtype;
  PetscTruth         opt;
  int                ierr;

  PetscFunctionBegin;
  /* Determine order of interpolation */
  if (numCorners == 3) {
    dtype = DISCRETIZATION_TRIANGULAR_2D_LINEAR;
  } else if (numCorners == 6) {
    dtype = DISCRETIZATION_TRIANGULAR_2D_QUADRATIC;
    ierr  = PetscOptionsHasName(mesh->prefix, "-mesh_move_linear", &opt);                                 CHKERRQ(ierr);
    if (opt == PETSC_TRUE) {
      dtype = DISCRETIZATION_TRIANGULAR_2D_LINEAR;
      /* Must play games here so that the coarse mesh does not get a MeshMover as well */
      ierr = PetscObjectReference((PetscObject) mover);                                                   CHKERRQ(ierr);
      ierr = MeshSetMover(mesh, PETSC_NULL);                                                              CHKERRQ(ierr);
      ierr = MeshCoarsen(mesh, PETSC_NULL, &movingMesh);                                                  CHKERRQ(ierr);
      ierr = MeshSetMover(mesh, mover);                                                                   CHKERRQ(ierr);
    }
  } else {
    SETERRQ1(PETSC_ERR_SUP, "Number of nodes per element %d not supported", numCorners);
  }

  /* Setup mesh velocity grid */
  ierr = GridCreate(movingMesh, &mover->meshVelocityGrid);                                                CHKERRQ(ierr);
  ierr = GridAddField(mover->meshVelocityGrid, "Velocity", dtype, dim, PETSC_NULL);                       CHKERRQ(ierr);
  ierr = GridAppendOptionsPrefix(mover->meshVelocityGrid, "mesh_vel_");                                   CHKERRQ(ierr);
  ierr = GridSetFromOptions(mover->meshVelocityGrid);                                                     CHKERRQ(ierr);
  ierr = GridSetActiveField(mover->meshVelocityGrid, 0);                                                  CHKERRQ(ierr);
  PetscLogObjectParent(movingMesh, mover->meshVelocityGrid);
  /* ALE operators need boundary values */
  mover->meshVelocityGrid->reduceElement = PETSC_TRUE;
  switch(mover->meshVelocityMethod) {
  case MESH_LAPLACIAN:
    ierr = GridSetMatOperator(mover->meshVelocityGrid, 0, 0, LAPLACIAN,    1.0, PETSC_FALSE, PETSC_NULL); CHKERRQ(ierr);
    break;
  case MESH_WEIGHTED_LAPLACIAN:
    ierr = GridSetMatOperator(mover->meshVelocityGrid, 0, 0, WEIGHTED_LAP, 1.0, PETSC_FALSE, PETSC_NULL); CHKERRQ(ierr);
    break;
  default:
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid velocity calculation method %d", mover->meshVelocityMethod);
  }

  /* Setup mesh acceleration grid */
  ierr = GridCreate(movingMesh, &mover->meshAccelerationGrid);                                            CHKERRQ(ierr);
  ierr = GridAddField(mover->meshAccelerationGrid, "Acceleration", dtype, dim, PETSC_NULL);               CHKERRQ(ierr);
  ierr = GridAppendOptionsPrefix(mover->meshAccelerationGrid, "mesh_accel_");                             CHKERRQ(ierr);
  ierr = GridSetFromOptions(mover->meshAccelerationGrid);                                                 CHKERRQ(ierr);
  ierr = GridSetActiveField(mover->meshAccelerationGrid, 0);                                              CHKERRQ(ierr);
  PetscLogObjectParent(mesh, mover->meshAccelerationGrid);
  switch(mover->meshAccelerationMethod) {
  case MESH_LAPLACIAN:
    ierr = GridSetMatOperator(mover->meshAccelerationGrid, 0, 0, LAPLACIAN, 1.0, PETSC_FALSE, PETSC_NULL);CHKERRQ(ierr);
    break;
  case MESH_WEIGHTED_LAPLACIAN:
    ierr = GridSetMatOperator(mover->meshVelocityGrid, 0, 0, WEIGHTED_LAP, 1.0, PETSC_FALSE, PETSC_NULL); CHKERRQ(ierr);
    break;
  default:
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid acceleration calculation method%d", mover->meshAccelerationMethod);
  }
  /* This will cause the coarse mesh to be destroyed when no longer needed */
  if (mesh != movingMesh) {
    ierr = MeshDestroy(movingMesh);                                                                       CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshMoverDestroy_Triangular_2D"
int MeshMoverDestroy_Triangular_2D(MeshMover mover) {
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshMapVertices_Triangular_2D"
int MeshMapVertices_Triangular_2D(GVec coarseVec, Mesh mesh, Vec vec) {
  Partition    p   = mesh->part;
  int          dim = mesh->dim;
  Grid         coarseGrid;
  Mesh         coarseMesh;
  Partition    coarseP;
  AO           coarseMap;
  PetscTruth   reduceSystem;
  int          numCoarseNodes;
  int         *classes;
  int          firstVar;
  int         *offsets, *localOffsets;
  int          reduceFirstVar     = 0;
  int         *reduceOffsets      = PETSC_NULL;
  int         *reduceLocalOffsets = PETSC_NULL;
  int        **reduceFieldClasses = PETSC_NULL;
  PetscScalar *vals, *coarseVals, *reduceVals;
  int          node, cNode, nclass, coarseVar, j;
  int          ierr;

  PetscFunctionBegin;
  ierr = GVecGetGrid(coarseVec,  &coarseGrid);                                                           CHKERRQ(ierr);
  ierr = GridGetMesh(coarseGrid, &coarseMesh);                                                           CHKERRQ(ierr);
  ierr = VecGetArray(vec,       &vals);                                                                  CHKERRQ(ierr);
  ierr = VecGetArray(coarseVec, &coarseVals);                                                            CHKERRQ(ierr);
  coarseP              = coarseMesh->part;
  coarseMap            = coarseMesh->coarseMap;
  reduceSystem         = coarseGrid->reduceSystem;
  numCoarseNodes       = coarseMesh->numNodes;
  reduceFieldClasses   = coarseGrid->cm->reduceFieldClasses;
  classes              = coarseGrid->cm->classes;
  firstVar             = coarseGrid->order->firstVar[p->rank];
  offsets              = coarseGrid->order->offsets;
  localOffsets         = coarseGrid->order->localOffsets;
  if (reduceSystem == PETSC_TRUE) {
    reduceOffsets      = coarseGrid->reduceOrder->offsets;
    reduceLocalOffsets = coarseGrid->reduceOrder->localOffsets;
    reduceFirstVar     = coarseGrid->reduceOrder->firstVar[p->rank];
    ierr = VecGetArray(coarseGrid->bdReduceVecCur, &reduceVals);                                         CHKERRQ(ierr);
  }
  for(cNode = 0; cNode < numCoarseNodes; cNode++) {
    /* Get corresponding node */
    ierr = PartitionLocalToGlobalNodeIndex(coarseP, cNode, &node);                                       CHKERRQ(ierr);
    ierr = AOPetscToApplication(coarseMap, 1, &node);                                                    CHKERRQ(ierr);
    ierr = PartitionGlobalToLocalNodeIndex(p,       node,  &node);                                       CHKERRQ(ierr);
    /* Calculate coarse offset */
    nclass = classes[cNode];
    if ((reduceSystem == PETSC_TRUE) && (reduceFieldClasses[0][nclass])) {
      if (cNode >= numCoarseNodes) {
        coarseVar = reduceLocalOffsets[cNode-numCoarseNodes];
      } else {
        coarseVar = reduceOffsets[cNode] - reduceFirstVar;
      }
      /* Set node value */
      for(j = 0; j < dim; j++) {
        vals[node*dim+j] = reduceVals[coarseVar+j];
      }
    } else {
      if (cNode >= numCoarseNodes) {
        coarseVar = localOffsets[cNode-numCoarseNodes];
      } else {
        coarseVar = offsets[cNode] - firstVar;
      }
      /* Set node value */
      for(j = 0; j < dim; j++) {
        vals[node*dim+j] = coarseVals[coarseVar+j];
      }
    }
  }
  ierr = VecRestoreArray(vec,       &vals);                                                              CHKERRQ(ierr);
  ierr = VecRestoreArray(coarseVec, &coarseVals);                                                        CHKERRQ(ierr);
  if (reduceSystem == PETSC_TRUE) {
    ierr = VecRestoreArray(coarseGrid->bdReduceVecCur, &reduceVals);                                     CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshInterpMidnodes_Triangular_2D"
int MeshInterpMidnodes_Triangular_2D(Mesh mesh, Vec vec) {
  Mesh_Triangular *tri         = (Mesh_Triangular *) mesh->data;
  int              dim         = mesh->dim;
  int              numCorners  = mesh->numCorners;
  int              numElements = mesh->numFaces;
  int             *elements    = tri->faces;
  PetscScalar     *vals;
  int              elem, corner, node, nNode1, nNode2, j;
  int              ierr;

  PetscFunctionBegin;
  ierr = VecGetArray(vec, &vals);                                                                        CHKERRQ(ierr);
  for(elem = 0; elem < numElements; elem++) {
    for(corner = 3; corner < numCorners; corner++) {
      node    = elements[elem*numCorners+corner];
      nNode1  = elements[elem*numCorners+((corner+1)%3)];
      nNode2  = elements[elem*numCorners+((corner+2)%3)];
      for(j = 0; j < dim; j++) {
        vals[node*dim+j] = 0.5*(vals[nNode1*dim+j] + vals[nNode2*dim+j]);
      }
    }
  }
  ierr = VecRestoreArray(vec, &vals);                                                                    CHKERRQ(ierr);
  PetscLogFlops(2*dim*(numCorners-3)*numElements);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshUpdateVelocities_Triangular_2D"
int MeshUpdateNodeValues_Triangular_2D(MeshMover mover, GVec sol, Vec vec, Vec ghostVec) {
  Mesh       mesh       = mover->mesh;
  Partition  p          = mesh->part;
  int        numCorners = mesh->numCorners;
  Grid       grid;
  int        numOverlapNodes, numFuncs;
  PetscTruth match;
  int        ierr;

  PetscFunctionBegin;
  /* Check grid */
  ierr = GVecGetGrid(sol, &grid);                                                                        CHKERRQ(ierr);
  ierr = PartitionGetNumOverlapNodes(p, &numOverlapNodes);                                               CHKERRQ(ierr);
  ierr = DiscretizationGetNumFunctions(grid->fields[0].disc, &numFuncs);                                 CHKERRQ(ierr);
  if (grid->reduceSystem == PETSC_TRUE) {
    if (numFuncs == numCorners) {
      if (grid->order->numOverlapVars + grid->reduceOrder->numOverlapVars != numOverlapNodes*2) {
        SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid number of mesh velocity unknowns %d should be %d",
                 grid->order->numOverlapVars + grid->reduceOrder->numOverlapVars, numOverlapNodes*2);
      }
    }
  } else {
    if (numFuncs == numCorners) {
      if (grid->order->numOverlapVars != numOverlapNodes*2) {
        SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid number of mesh velocity unknowns %d should be %d",
                 grid->order->numOverlapVars, numOverlapNodes*2);
      }
    }
  }
  /* Map coarse nodes to fine vertices */
  ierr = PetscTypeCompare((PetscObject) grid->fields[0].disc, DISCRETIZATION_TRIANGULAR_2D_LINEAR, &match);CHKERRQ(ierr);
  if ((numCorners == 6) && (match == PETSC_TRUE)) {
    ierr = MeshMapVertices_Triangular_2D(sol, mesh, vec);                                                 CHKERRQ(ierr);
  } else {
    ierr = VecCopy(sol, vec);                                                                             CHKERRQ(ierr);
  }
  /* Copy to ghost vector */
  ierr = VecScatterBegin(vec, ghostVec, INSERT_VALUES, SCATTER_FORWARD, mover->meshVelocityScatter);      CHKERRQ(ierr);
  ierr = VecScatterEnd(vec, ghostVec, INSERT_VALUES, SCATTER_FORWARD, mover->meshVelocityScatter);        CHKERRQ(ierr);
  /* Interpolate velocities of midnodes */
  ierr = MeshInterpMidnodes_Triangular_2D(mesh, ghostVec);                                                CHKERRQ(ierr);

  /* REMOVE THIS AFTER YOU REWRITE THE ALE STUFF */
  if (sol == mover->meshVelocitySol) {
    ierr = GridGlobalToLocal(mover->meshVelocityGrid, INSERT_VALUES, sol);                                CHKERRQ(ierr);
  } else if (sol == mover->meshAccelerationSol) {
    ierr = GridGlobalToLocal(mover->meshAccelerationGrid, INSERT_VALUES, sol);                            CHKERRQ(ierr);
  } else {
    SETERRQ(PETSC_ERR_SUP, "Unrecognized solution vector");
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshUpdateCoords_Triangular_2D"
int MeshUpdateCoords_Triangular_2D(Mesh mesh, PetscScalar *v, PetscScalar *a, double t) {
  Mesh_Triangular *tri   = (Mesh_Triangular *) mesh->data;
  Partition        p     = mesh->part;
  int              dim   = mesh->dim;
  double          *nodes = tri->nodes;
  int              numOverlapNodes;
  int              node;
  int              ierr;

  PetscFunctionBegin;
  ierr = PartitionGetNumOverlapNodes(p, &numOverlapNodes);                                               CHKERRQ(ierr);
  /* x' = x + v t + 0.5 a t^2 */
  if (mesh->isPeriodic == PETSC_TRUE) {
    for(node = 0; node < numOverlapNodes; node++) {
      nodes[node*2]   = MeshPeriodicX(mesh, nodes[node*2]   + PetscRealPart(v[node*2])*t   + 0.5*PetscRealPart(a[node*2])*t*t);
      nodes[node*2+1] = MeshPeriodicY(mesh, nodes[node*2+1] + PetscRealPart(v[node*2+1])*t + 0.5*PetscRealPart(a[node*2+1])*t*t);
#ifdef DEBUG_BD_MOVEMENT
      if (tri->markers[node] > 0)
        PetscPrintf(PETSC_COMM_SELF, "node %d: u = %lf v = %lf\n", node, v[node*2], v[node*2+1]);
#endif
    }
  } else {
    for(node = 0; node < numOverlapNodes*dim; node++) {
      nodes[node] += PetscRealPart(v[node])*t + 0.5*PetscRealPart(a[node])*t*t;
#ifdef DEBUG_BD_MOVEMENT
      if ((node%2 == 0) && (tri->markers[node/2] > 0))
        PetscPrintf(PETSC_COMM_SELF, "node %d: u = %lf v = %lf\n", node/2, v[node], v[node+1]);
#endif
    }
  }
  PetscLogFlops(6*numOverlapNodes*dim);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSyncCoords_Triangular_2D"
int MeshSyncCoords_Triangular_2D(Mesh mesh, Mesh coarseMesh) {
  Mesh_Triangular *tri            = (Mesh_Triangular *) mesh->data;
  Mesh_Triangular *coarseTri      = (Mesh_Triangular *) coarseMesh->data;
  Partition        p              = mesh->part;
  Partition        coarseP        = coarseMesh->part;
  AO               coarseMap      = coarseMesh->coarseMap;
  int              dim            = mesh->dim;
  int              numCoarseNodes = coarseMesh->numNodes;
  double          *nodes          = tri->nodes;
  double          *coarseNodes    = coarseTri->nodes;
  int              node, cNode, j;
  int              ierr;

  PetscFunctionBegin;
  for(cNode = 0; cNode < numCoarseNodes; cNode++) {
    ierr = PartitionLocalToGlobalNodeIndex(coarseP, cNode, &node);                                       CHKERRQ(ierr);
    ierr = AOPetscToApplication(coarseMap, 1, &node);                                                    CHKERRQ(ierr);
    ierr = PartitionGlobalToLocalNodeIndex(p,       node,  &node);                                       CHKERRQ(ierr);
    for(j = 0; j < dim; j++) {
      coarseNodes[cNode*dim+j] = nodes[node*dim+j];
    }
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshMoveMesh_Triangular_2D"
int MeshMoveMesh_Triangular_2D(MeshMover mover, double t) {
  Mesh         mesh    = mover->mesh;
  Mesh         velMesh = mover->meshVelocityGrid->mesh;
  PetscScalar *v, *a;
  PetscDraw    draw;
  int          ierr;

  PetscFunctionBegin;
  ierr = VecGetArray(mover->nodeVelocitiesGhost,    &v);                                                  CHKERRQ(ierr);
  ierr = VecGetArray(mover->nodeAccelerationsGhost, &a);                                                  CHKERRQ(ierr);

  /* Move mesh nodes */
  ierr = MeshUpdateCoords_Triangular_2D(mesh, v, a, t);                                                   CHKERRQ(ierr);
  /* Move mesh velocity and acceleration mesh nodes */
  ierr = MeshSyncCoords_Triangular_2D(mesh, velMesh);                                                     CHKERRQ(ierr);

  ierr = VecRestoreArray(mover->nodeVelocitiesGhost,    &v);                                              CHKERRQ(ierr);
  ierr = VecRestoreArray(mover->nodeAccelerationsGhost, &a);                                              CHKERRQ(ierr);
 
  /* Update bounding rectangle */
  ierr = MeshUpdateBoundingBox(mesh);                                                                     CHKERRQ(ierr);

  /* View mesh */
  if (mover->movingMeshViewer) {
    ierr = PetscViewerDrawGetDraw(mover->movingMeshViewer, 0, &draw);                                     CHKERRQ(ierr);
    ierr = PetscDrawClear(draw);                                                                          CHKERRQ(ierr);
    ierr = PetscDrawSetTitle(draw, mover->movingMeshViewerCaption);                                       CHKERRQ(ierr);
    ierr = MeshView(mesh, mover->movingMeshViewer);                                                       CHKERRQ(ierr);
    ierr = PetscViewerFlush(mover->movingMeshViewer);                                                     CHKERRQ(ierr);
    ierr = PetscDrawPause(draw);                                                                          CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshReform_Triangular_2D"
/*
  MeshReform_Triangular_2D - Reforms a two dimensional unstructured mesh
  using the original boundary.

  Input Parameters:
. mesh    - The mesh begin reformed
. refine  - Flag indicating whether to refine or recreate the mesh
. area    - A function which gives an area constraint when evaluated inside an element
. newBd   - Flag to determine whether the boundary should be generated (PETSC_TRUE) or taken from storage

  Output Parameter:
. newMesh - The reformed mesh

.keywords unstructured grid
.seealso GridReform(), MeshRefine(), MeshSetBoundary()
*/
int MeshReform_Triangular_2D(Mesh mesh, PetscTruth refine, PointFunction area, PetscTruth newBd, Mesh *newMesh) {
  Mesh                     m;
  Mesh_Triangular         *tri = (Mesh_Triangular *) mesh->data;
  Partition                p   = mesh->part;
  Partition_Triangular_2D *q   = (Partition_Triangular_2D *) p->data;
  MeshBoundary2D           bdCtx;
  int                      numEdges     = mesh->numBdEdges;
  int                      origLocNodes = mesh->numBdNodes;
  int                      numLocNodes;
  int                     *firstBdEdge;
  int                     *firstBdNode;
  int                     *firstHole;
  int                     *numRecvNodes;
  int                     *numRecvMarkers;
  int                     *numRecvSegments;
  int                     *numRecvSegMarkers;
  int                     *numRecvHoles;
  int                     *nodeOffsets;
  int                     *markerOffsets;
  int                     *segmentOffsets;
  int                     *segMarkerOffsets;
  int                     *holeOffsets;
  double                  *locNodes;
  int                     *locMarkers;
  int                     *locSegments;
  int                     *locSegMarkers;
  double                  *nodes;
  int                     *markers;
  int                     *segments;
  int                     *segMarkers;
  double                  *holes;
  int                     *globalNodes;
  int                     *aOrder;
  int                      rank, numProcs;
  int                      numLocSegments, numSegments, numNodes, numHoles;
  int                      proc, bd, gNode, lNode, bdNode, edge, bdEdge, edgeEnd;
  int                      ierr;

  PetscFunctionBegin;
  ierr = MPI_Comm_size(mesh->comm, &numProcs);                                                            CHKERRQ(ierr);
  ierr = MPI_Comm_rank(mesh->comm, &rank);                                                                CHKERRQ(ierr);
  if ((newBd == PETSC_TRUE) || (mesh->bdCtxNew == PETSC_NULL)) {
    ierr = MPI_Reduce(&mesh->numHoles, &numHoles, 1, MPI_INT, MPI_SUM, 0, mesh->comm);                    CHKERRQ(ierr);
    ierr = PetscMalloc((numProcs+1)   * sizeof(int),    &firstBdEdge);                                    CHKERRQ(ierr);
    ierr = PetscMalloc((numProcs+1)   * sizeof(int),    &firstBdNode);                                    CHKERRQ(ierr);
    ierr = PetscMalloc((numProcs+1)   * sizeof(int),    &firstHole);                                      CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &numRecvNodes);                                   CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &numRecvMarkers);                                 CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &numRecvSegments);                                CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &numRecvSegMarkers);                              CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &numRecvHoles);                                   CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &nodeOffsets);                                    CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &markerOffsets);                                  CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &segmentOffsets);                                 CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &segMarkerOffsets);                               CHKERRQ(ierr);
    ierr = PetscMalloc(numProcs       * sizeof(int),    &holeOffsets);                                    CHKERRQ(ierr);
    ierr = PetscMalloc(origLocNodes*2 * sizeof(double), &locNodes);                                       CHKERRQ(ierr);
    ierr = PetscMalloc(origLocNodes   * sizeof(int),    &locMarkers);                                     CHKERRQ(ierr);
    ierr = PetscMalloc(numEdges*2     * sizeof(int),    &locSegments);                                    CHKERRQ(ierr);
    ierr = PetscMalloc(numEdges       * sizeof(int),    &locSegMarkers);                                  CHKERRQ(ierr);
    ierr = PetscMalloc(origLocNodes   * sizeof(int),    &globalNodes);                                    CHKERRQ(ierr);

    /* Make all node numbers global */
    numLocNodes    = 0;
    numLocSegments = 0;
    for(bd = 0; bd < mesh->numBd; bd++) {
      /* Get boundary nodes */
      for(bdNode = tri->bdBegin[bd]; bdNode < tri->bdBegin[bd+1]; bdNode++) {
        gNode = tri->bdNodes[bdNode];
        lNode = gNode - q->firstNode[p->rank];
        if ((lNode >= 0) && (lNode < mesh->numNodes)) {
          if (tri->degrees[lNode] > 0) {
            locNodes[numLocNodes*2]   = tri->nodes[lNode*2];
            locNodes[numLocNodes*2+1] = tri->nodes[lNode*2+1];
            locMarkers[numLocNodes]   = tri->bdMarkers[bd];
            globalNodes[numLocNodes]  = gNode;
            numLocNodes++;
          }
        }
      }

      /* Get boundary segments */
      for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) {
        edge = tri->bdEdges[bdEdge] - q->firstEdge[rank];
        if ((edge >= 0) && (edge < mesh->numEdges)) {
          for(edgeEnd = 0; edgeEnd < 2; edgeEnd++) {
            lNode = tri->edges[edge*2+edgeEnd];
            ierr  = PartitionLocalToGlobalNodeIndex(p, lNode, &gNode);                                    CHKERRQ(ierr);
            locSegments[numLocSegments*2+edgeEnd] = gNode;
          }
          locSegMarkers[numLocSegments]   = tri->bdMarkers[bd];
          numLocSegments++;
        }
      }
    }

    /* Assemble offsets */
    ierr = MPI_Allgather(&numLocSegments, 1, MPI_INT, &firstBdEdge[1], 1, MPI_INT, p->comm);              CHKERRQ(ierr);
    firstBdEdge[0] = 0;
    for(proc = 1; proc <= numProcs; proc++)
      firstBdEdge[proc] = firstBdEdge[proc] + firstBdEdge[proc-1];
    numSegments = firstBdEdge[numProcs];
    if ((rank == 0) && (numSegments != numEdges)) SETERRQ(PETSC_ERR_PLIB, "Invalid number of boundary segments");
    ierr = MPI_Allgather(&numLocNodes,    1, MPI_INT, &firstBdNode[1], 1, MPI_INT, p->comm);              CHKERRQ(ierr);
    firstBdNode[0] = 0;
    for(proc = 1; proc <= numProcs; proc++)
      firstBdNode[proc] = firstBdNode[proc] + firstBdNode[proc-1];
    numNodes = firstBdNode[numProcs];
    ierr = MPI_Allgather(&mesh->numHoles,  1, MPI_INT, &firstHole[1],   1, MPI_INT, p->comm);             CHKERRQ(ierr);
    firstHole[0] = 0;
    for(proc = 1; proc <= numProcs; proc++)
      firstHole[proc] = firstHole[proc] + firstHole[proc-1];
    if ((rank == 0) && (numHoles != firstHole[numProcs])) SETERRQ(PETSC_ERR_PLIB, "Invalid number of holes");

    for(proc = 0; proc < numProcs; proc++) {
      numRecvNodes[proc]      = (firstBdNode[proc+1] - firstBdNode[proc])*2;
      numRecvMarkers[proc]    = (firstBdNode[proc+1] - firstBdNode[proc]);
      numRecvSegments[proc]   = (firstBdEdge[proc+1] - firstBdEdge[proc])*2;
      numRecvSegMarkers[proc] = (firstBdEdge[proc+1] - firstBdEdge[proc]);
      numRecvHoles[proc]      = (firstHole[proc+1]   - firstHole[proc])*2;
    }
    nodeOffsets[0]      = 0;
    markerOffsets[0]    = 0;
    segmentOffsets[0]   = 0;
    segMarkerOffsets[0] = 0;
    holeOffsets[0]      = 0;
    for(proc = 1; proc < numProcs; proc++) {
      nodeOffsets[proc]      = numRecvNodes[proc-1]      + nodeOffsets[proc-1];
      markerOffsets[proc]    = numRecvMarkers[proc-1]    + markerOffsets[proc-1];
      segmentOffsets[proc]   = numRecvSegments[proc-1]   + segmentOffsets[proc-1];
      segMarkerOffsets[proc] = numRecvSegMarkers[proc-1] + segMarkerOffsets[proc-1];
      holeOffsets[proc]      = numRecvHoles[proc-1]      + holeOffsets[proc-1];
    }

    if (rank == 0) {
      ierr = PetscMalloc(numNodes*2 * sizeof(double), &nodes);                                            CHKERRQ(ierr);
      ierr = PetscMalloc(numNodes   * sizeof(int),    &markers);                                          CHKERRQ(ierr);
      ierr = PetscMalloc(numEdges*2 * sizeof(int),    &segments);                                         CHKERRQ(ierr);
      ierr = PetscMalloc(numEdges   * sizeof(int),    &segMarkers);                                       CHKERRQ(ierr);
      if (numHoles) {
        ierr  = PetscMalloc(numHoles*2 * sizeof(double), &holes);                                         CHKERRQ(ierr);
      } else {
        holes = PETSC_NULL;
      }
    }

    /* Make new numbering for boundary nodes */
    ierr = PetscMalloc(numNodes * sizeof(int), &aOrder);                                                  CHKERRQ(ierr);
    ierr = MPI_Allgatherv(globalNodes, numLocNodes, MPI_INT, aOrder, numRecvMarkers, markerOffsets, MPI_INT, p->comm);
                                                                                                          CHKERRQ(ierr);

    /* Renumber segments */
    for(edge = 0; edge < numLocSegments*2; edge++) {
      gNode = locSegments[edge];
      for(bdNode = 0; bdNode < numNodes; bdNode++)
        if (aOrder[bdNode] == gNode)
          break;
      if (bdNode == numNodes) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary node");
      locSegments[edge] = bdNode;
    }
    ierr = PetscFree(aOrder);                                                                             CHKERRQ(ierr);

    ierr = MPI_Gatherv(locNodes,      numLocNodes*2,    MPI_DOUBLE, nodes,      numRecvNodes,      nodeOffsets,      MPI_DOUBLE, 0, mesh->comm);
                                                                                                          CHKERRQ(ierr);
    ierr = MPI_Gatherv(locMarkers,    numLocNodes,      MPI_INT,    markers,    numRecvMarkers,    markerOffsets,    MPI_INT,    0, mesh->comm);
                                                                                                          CHKERRQ(ierr);
    ierr = MPI_Gatherv(locSegments,   numLocSegments*2, MPI_INT,    segments,   numRecvSegments,   segmentOffsets,   MPI_INT,    0, mesh->comm);
                                                                                                          CHKERRQ(ierr);
    ierr = MPI_Gatherv(locSegMarkers, numLocSegments,   MPI_INT,    segMarkers, numRecvSegMarkers, segMarkerOffsets, MPI_INT,    0, mesh->comm);
                                                                                                          CHKERRQ(ierr);
    ierr = MPI_Gatherv(mesh->holes,   mesh->numHoles*2, MPI_DOUBLE, holes,      numRecvHoles,      holeOffsets,      MPI_DOUBLE, 0, mesh->comm);
                                                                                                          CHKERRQ(ierr);

    /* Cleanup */
    ierr = PetscFree(firstBdEdge);                                                                        CHKERRQ(ierr);
    ierr = PetscFree(firstBdNode);                                                                        CHKERRQ(ierr);
    ierr = PetscFree(firstHole);                                                                          CHKERRQ(ierr);
    ierr = PetscFree(numRecvNodes);                                                                       CHKERRQ(ierr);
    ierr = PetscFree(numRecvMarkers);                                                                     CHKERRQ(ierr);
    ierr = PetscFree(numRecvSegments);                                                                    CHKERRQ(ierr);
    ierr = PetscFree(numRecvSegMarkers);                                                                  CHKERRQ(ierr);
    ierr = PetscFree(numRecvHoles);                                                                       CHKERRQ(ierr);
    ierr = PetscFree(nodeOffsets);                                                                        CHKERRQ(ierr);
    ierr = PetscFree(markerOffsets);                                                                      CHKERRQ(ierr);
    ierr = PetscFree(segmentOffsets);                                                                     CHKERRQ(ierr);
    ierr = PetscFree(segMarkerOffsets);                                                                   CHKERRQ(ierr);
    ierr = PetscFree(holeOffsets);                                                                        CHKERRQ(ierr);
    ierr = PetscFree(locNodes);                                                                           CHKERRQ(ierr);
    ierr = PetscFree(locMarkers);                                                                         CHKERRQ(ierr);
    ierr = PetscFree(locSegments);                                                                        CHKERRQ(ierr);
    ierr = PetscFree(locSegMarkers);                                                                      CHKERRQ(ierr);
    ierr = PetscFree(globalNodes);                                                                        CHKERRQ(ierr);

    bdCtx.numVertices = numNodes;
    bdCtx.vertices    = nodes;
    bdCtx.markers     = markers;
    bdCtx.numSegments = numEdges;
    bdCtx.segments    = segments;
    bdCtx.segMarkers  = segMarkers;
    bdCtx.numHoles    = numHoles;
    bdCtx.holes       = holes;
  } else {
    bdCtx.numVertices = mesh->bdCtxNew->numVertices;
    bdCtx.vertices    = mesh->bdCtxNew->vertices;
    bdCtx.markers     = mesh->bdCtxNew->markers;
    bdCtx.numSegments = mesh->bdCtxNew->numSegments;
    bdCtx.segments    = mesh->bdCtxNew->segments;
    bdCtx.segMarkers  = mesh->bdCtxNew->segMarkers;
    bdCtx.numHoles    = mesh->bdCtxNew->numHoles;
    bdCtx.holes       = mesh->bdCtxNew->holes;
  }
  bdCtx.numBd = mesh->numBd;

  ierr = MeshCreate(mesh->comm, &m);                                                                      CHKERRQ(ierr);
  ierr = MeshSetDimension(m, 2);                                                                          CHKERRQ(ierr);
  ierr = MeshSetPeriodicDimension(m, 0, mesh->isPeriodicDim[0]);                                          CHKERRQ(ierr);
  ierr = MeshSetPeriodicDimension(m, 1, mesh->isPeriodicDim[1]);                                          CHKERRQ(ierr);
  ierr = MeshSetNumCorners(m, mesh->numCorners);                                                          CHKERRQ(ierr);
  ierr = MeshSetBoundary(m, &bdCtx);                                                                      CHKERRQ(ierr);
  ierr = MeshSetType(m, mesh->type_name);                                                                 CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) m, "Mesh");                                                     CHKERRQ(ierr);
  *newMesh = m;

  if (rank == 0) {
    if ((newBd == PETSC_TRUE) || (mesh->bdCtxNew == PETSC_NULL)) {
      ierr = PetscFree(nodes);                                                                            CHKERRQ(ierr);
      ierr = PetscFree(markers);                                                                          CHKERRQ(ierr);
      ierr = PetscFree(segments);                                                                         CHKERRQ(ierr);
      ierr = PetscFree(segMarkers);                                                                       CHKERRQ(ierr);
      if (numHoles) {
        ierr = PetscFree(holes);                                                                          CHKERRQ(ierr);
      }
    } else {
#if 0
      ierr = MeshBoundaryDestroy(mesh->bdCtxNew);                                                         CHKERRQ(ierr);
#endif
      mesh->bdCtxNew = PETSC_NULL;
    }
  }

#ifdef PETSC_USE_BOPT_g
  ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                           CHKERRQ(ierr);
#endif
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshDebugBoundary_Triangular_2D_Private"
int MeshDebugBoundary_Triangular_2D_Private(Mesh mesh, GVec nodeVelocities) {
  Mesh_Triangular *tri      = (Mesh_Triangular *) mesh->data;
  int              numNodes = mesh->numNodes;
  double          *nodes    = tri->nodes;
  int             *markers  = tri->markers;
  int              rank     = mesh->part->rank;
  PetscScalar     *vel;
  int              node;
  int              ierr;

  PetscFunctionBegin;
  ierr = VecGetArray(nodeVelocities, &vel);                                                              CHKERRQ(ierr);
  for(node = 0; node < numNodes; node++) {
    if (markers[node] > 0) {
      PetscPrintf(PETSC_COMM_SELF, "[%d]node %d at (%g,%g) vel: (%g,%g)\n",
                  rank, node, nodes[node*2], nodes[node*2+1], PetscRealPart(vel[node*2]), PetscRealPart(vel[node*2+1]));
    }
  }
  ierr = VecRestoreArray(nodeVelocities, &vel);                                                          CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshDebugBoundary_Triangular_2D"
int MeshDebugBoundary_Triangular_2D(MeshMover mover) {
  Mesh             mesh               = mover->mesh;
  Mesh_Triangular *tri                = (Mesh_Triangular *) mesh->data;
  Grid             velGrid            = mover->meshVelocityGrid;
  Grid             accGrid            = mover->meshAccelerationGrid;
  int              numBd              = mesh->numBd;
  int             *bdMarkers          = tri->bdMarkers;
  int              numNodes           = mesh->numNodes;
  int              rank               = mesh->part->rank;
  int            **reduceFieldClasses = velGrid->cm->reduceFieldClasses;
  int             *classes            = velGrid->cm->classes;
  int              firstVar           = velGrid->order->firstVar[rank];
  int             *offsets            = velGrid->order->offsets;
  int             *localOffsets       = velGrid->order->localOffsets;
  int              reduceFirstVar     = 0;
  int             *reduceOffsets      = PETSC_NULL;
  int             *reduceLocalOffsets = PETSC_NULL;
  PetscScalar     *v, *reduceV, *reduceA;
  int              b, bd, node, nclass, var;
  int              ierr;

  PetscFunctionBegin;
  ierr = VecGetArray(velGrid->ghostVec, &v);                                                             CHKERRQ(ierr);
#ifdef DEBUG_ACC
  ierr = VecGetArray(accGrid->ghostVec, &a);                                                             CHKERRQ(ierr);
#endif
  if (velGrid->reduceSystem == PETSC_TRUE) {
    reduceOffsets      = velGrid->reduceOrder->offsets;
    reduceLocalOffsets = velGrid->reduceOrder->localOffsets;
    reduceFirstVar     = velGrid->reduceOrder->firstVar[rank];
    ierr = VecGetArray(velGrid->bdReduceVecCur, &reduceV);                                               CHKERRQ(ierr);
    ierr = VecGetArray(accGrid->bdReduceVecCur, &reduceA);                                               CHKERRQ(ierr);
  }
  for(b = 0; b < numBd; b++) {
    bd = bdMarkers[b];
    if (bd < 0) continue;

    ierr = MeshGetBoundaryStart(mesh, bd, PETSC_TRUE, &node);                                            CHKERRQ(ierr);
    while (node >= 0) {
      /* Print the value of velocity on boundary node */
      nclass = classes[node];
      if ((reduceFieldClasses != PETSC_NULL) && (reduceFieldClasses[0][nclass])) {
        if (node >= numNodes) {
          var = reduceLocalOffsets[node-numNodes];
        } else {
          var = reduceOffsets[node] - reduceFirstVar;
        }
        PetscPrintf(PETSC_COMM_SELF, "[%d]node %d: u = %g v = %g", rank, node, PetscRealPart(reduceV[var]), PetscRealPart(reduceV[var+1]));
#ifdef DEBUG_ACC
        PetscPrintf(PETSC_COMM_SELF, "[%d]node %d: a = %g b = %g", rank, node, PetscRealPart(reduceA[var]), PetscRealPart(reduceA[var+1]));
#endif
      } else {
        if (node >= numNodes) {
          var = localOffsets[node-numNodes];
        } else {
          var = offsets[node] - firstVar;
        }
        PetscPrintf(PETSC_COMM_SELF, "[%d]node %d: u = %g v = %g", rank, node, PetscRealPart(v[var]), PetscRealPart(v[var+1]));
#ifdef DEBUG_ACC
        PetscPrintf(PETSC_COMM_SELF, "[%d]node %d: a = %g b = %g", rank, node, PetscRealPart(a[var]), PetscRealPart(a[var+1]));
#endif
      }
      ierr = MeshGetBoundaryNext(mesh, bd, PETSC_TRUE, &node);                                           CHKERRQ(ierr);
    }
  }
  ierr = VecRestoreArray(velGrid->ghostVec, &v);                                                         CHKERRQ(ierr);
#ifdef DEBUG_ACC
  ierr = VecRestoreArray(accGrid->ghostVec, &a);                                                         CHKERRQ(ierr);
#endif
  if (velGrid->reduceSystem == PETSC_TRUE) {
    ierr = VecRestoreArray(velGrid->bdReduceVecCur, &reduceV);                                           CHKERRQ(ierr);
    ierr = VecRestoreArray(accGrid->bdReduceVecCur, &reduceA);                                           CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

static struct _MeshMoverOps MeshMoverOps = {/* Generic Operations */
                                            MeshMoverSetup_Triangular_2D,
                                            PETSC_NULL/* MeshMoverSetFromOptions */,
                                            MeshMoverDestroy_Triangular_2D,
                                            /* Mesh-Specific Operations */
                                            MeshMoveMesh_Triangular_2D,
                                            MeshUpdateNodeValues_Triangular_2D,
                                            MeshReform_Triangular_2D};

EXTERN_C_BEGIN
#undef  __FUNCT__
#define __FUNCT__ "MeshMoverCreate_Triangular_2D"
int MeshMoverCreate_Triangular_2D(MeshMover mover, ParameterDict dict) {
  int ierr;

  PetscFunctionBegin;
  ierr = PetscMemcpy(mover->ops, &MeshMoverOps, sizeof(struct _MeshMoverOps));                            CHKERRQ(ierr);

  PetscFunctionReturn(0);
}
EXTERN_C_END