#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: meshQuery.c,v 1.19 2000/10/17 13:48:57 knepley Exp $";
#endif

/*
     Defines the interface to the mesh query functions
*/

#include "src/mesh/meshimpl.h"    /*I "mesh.h" I*/

/*----------------------------------------------- Mesh Query Functions -----------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "MeshSetDimension"
/*@
  MeshSetDimension - This function sets the dimension of the Mesh. This may only be done before a mesh is generated.

  Not collective

  Input Parameters:
+ mesh - The mesh
- dim  - The mesh dimension

  Level: intermediate

.keywords: Mesh, dimension
.seealso: MatGetDimension(), PartitionGetDimension()
@*/
int MeshSetDimension(Mesh mesh, int dim)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if ((dim < 0) || (dim > 3)) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Dimension %d is invalid", dim);
  if (mesh->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Dimension cannot be set after mesh has been generated");
  mesh->dim = dim;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetDimension"
/*@
  MeshGetDimension - This function returns the dimension of the Mesh

  Not collective

  Input Parameter:
. mesh - The mesh

  Output Parameter:
. dim  - The mesh dimension

  Level: intermediate

.keywords: Mesh, dimension
.seealso: MatSetDimension(), PartitionGetDimension()
@*/
int MeshGetDimension(Mesh mesh, int *dim)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(dim);
  *dim = mesh->dim;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetInfo"
/*@
  MeshGetInfo - This function returns some information about the mesh.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Output Parameters:
+ vertices - The number of vertices
. nodes    - The number of nodes
. edges    - The number of edges
- elements - The number of elements

  Level: intermediate

.keywords: mesh
.seealso: MatGetCorners()
@*/
int MeshGetInfo(Mesh mesh, int *vertices, int *nodes, int *edges, int *elements) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (vertices != PETSC_NULL) *vertices = mesh->numVertices;
  if (nodes    != PETSC_NULL) *nodes    = mesh->numNodes;
  if (edges    != PETSC_NULL) *edges    = mesh->numEdges;
  if (elements != PETSC_NULL) *elements = mesh->numFaces;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSetNumCorners"
/*@
  MeshSetNumCorners - This function sets the number of nodes on each element. This may only be done before a mesh is generated.

  Not collective

  Input Parameters:
+ mesh       - The mesh
- numCorners - The number of nodes on each element

  Level: intermediate

.keywords: Mesh, dimension
.seealso: MatGetNumCorners()
@*/
int MeshSetNumCorners(Mesh mesh, int numCorners)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (numCorners <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Number of corners %d is invalid", numCorners);
  if (mesh->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Corners cannot be set after mesh has been generated");
  mesh->numCorners = numCorners;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNumCorners"
/*@
  MeshGetNumCorners - This function returns the number of nodes on each element.

  Not collective

  Input Parameter:
. mesh       - The mesh

  Output Parameter:
. numCorners - The number of nodes on an element

  Level: intermediate

.keywords: mesh, corner
.seealso: MatGetInfo()
@*/
int MeshGetNumCorners(Mesh mesh, int *numCorners) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(numCorners);
  *numCorners = mesh->numCorners;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSetBoundingBox"
/*@
  MeshSetBoundingBox - This function sets the bounding box for the mesh. This can only be done before
  the mesh is generated.

  Not collective

  Input Parameters:
+ mesh                   - The mesh
. startX, startY, startZ - The lower-left corner of the box
- endX,   endY,   endZ   - The upper-right corner of the box

  Level: intermediate

.keywords: mesh, bounding box
.seealso: MeshGetBoundingBox(), MeshSetLocalBoundingBox(), MeshUpdateBoundingBox()
@*/
int MeshSetBoundingBox(Mesh mesh, PetscReal startX, PetscReal startY, PetscReal startZ, PetscReal endX, PetscReal endY, PetscReal endZ)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if ((startX > endX) || (startY > endY) || (startZ > endZ)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "The bounding box is inverted");
  if (mesh->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Bounding box cannot be set after mesh has been generated");
  mesh->startX = startX;
  mesh->startY = startY;
  mesh->startZ = startZ;
  mesh->endX   = endX;
  mesh->endY   = endY;
  mesh->endZ   = endZ;
  mesh->sizeX  = endX - startX;
  mesh->sizeY  = endY - startY;
  mesh->sizeZ  = endZ - startZ;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetBoundingBox"
/*@
  MeshGetBoundingBox - This function returns the bounding box for the mesh

  Not collective

  Input Parameter:
. mesh - The mesh

  Output Parameters:
+ startX, startY, startZ - The lower-left corner of the box
- endX,   endY,   endZ   - The upper-right corner of the box

  Level: intermediate

.keywords: mesh, bounding box
.seealso: MeshSetBoundingBox(), MeshGetLocalBoundingBox(), MeshUpdateBoundingBox()
@*/
int MeshGetBoundingBox(Mesh mesh, PetscReal *startX, PetscReal *startY, PetscReal *startZ, PetscReal *endX, PetscReal *endY, PetscReal *endZ)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (startX != PETSC_NULL) {
    PetscValidPointer(startX);
    *startX = mesh->startX;
  }
  if (startY != PETSC_NULL) {
    PetscValidPointer(startY);
    *startY = mesh->startY;
  }
  if (startZ != PETSC_NULL) {
    PetscValidPointer(startZ);
    *startZ = mesh->startZ;
  }
  if (endX   != PETSC_NULL) {
    PetscValidPointer(endX);
    *endX   = mesh->endX;
  }
  if (endY   != PETSC_NULL) {
    PetscValidPointer(endY);
    *endY   = mesh->endY;
  }
  if (endZ   != PETSC_NULL) {
    PetscValidPointer(endZ);
    *endZ   = mesh->endZ;
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSetLocalBoundingBox"
/*@
  MeshSetLocalBoundingBox - This function sets the local bounding box for the mesh. This can only be done before
  the mesh is generated. The local box is the smallest rectangle enclosing the portion of the mesh allocated to
  this processor.

  Not collective

  Input Parameters:
+ mesh                   - The mesh
. startX, startY, startZ - The lower-left corner of the local box
- endX,   endY,   endZ   - The upper-right corner of the local box

  Level: intermediate

.keywords: mesh, bounding box, local
.seealso: MeshGetLocalBoundingBox(), MeshSetBoundingBox(), MeshUpdateBoundingBox()
@*/
int MeshSetLocalBoundingBox(Mesh mesh, PetscReal startX, PetscReal startY, PetscReal startZ, PetscReal endX, PetscReal endY, PetscReal endZ)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if ((startX > endX) || (startY > endY) || (startZ > endZ)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "The bounding box is inverted");
  if (mesh->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Bounding box cannot be set after mesh has been generated");
  mesh->locStartX = startX;
  mesh->locStartY = startY;
  mesh->locStartZ = startZ;
  mesh->locEndX   = endX;
  mesh->locEndY   = endY;
  mesh->locEndZ   = endZ;
  mesh->locSizeX  = mesh->locEndX - mesh->locStartX;
  mesh->locSizeY  = mesh->locEndY - mesh->locStartY;
  mesh->locSizeZ  = mesh->locEndZ - mesh->locStartZ;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetLocalBoundingBox"
/*@
  MeshGetLocalBoundingBox - This function returns the local bounding box for the mesh. The local box is the smallest
  rectangle enclosing the portion of the mesh allocated to this processor.

  Not collective

  Input Parameter:
. mesh - The mesh

  Output Parameters:
+ startX, startY, startZ - The lower-left corner of the local box
- endX,   endY,   endZ   - The upper-right corner of the local box

  Level: intermediate

.keywords: mesh, bounding box, local
.seealso: MeshSetLocalBoundingBox(), MeshGetBoundingBox(), MeshUpdateBoundingBox()
@*/
int MeshGetLocalBoundingBox(Mesh mesh, PetscReal *startX, PetscReal *startY, PetscReal *startZ, PetscReal *endX, PetscReal *endY, PetscReal *endZ)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (startX != PETSC_NULL) {
    PetscValidPointer(startX);
    *startX = mesh->locStartX;
  }
  if (startY != PETSC_NULL) {
    PetscValidPointer(startY);
    *startY = mesh->locStartY;
  }
  if (startZ != PETSC_NULL) {
    PetscValidPointer(startZ);
    *startZ = mesh->locStartZ;
  }
  if (endX   != PETSC_NULL) {
    PetscValidPointer(endX);
    *endX   = mesh->locEndX;
  }
  if (endY   != PETSC_NULL) {
    PetscValidPointer(endY);
    *endY   = mesh->locEndY;
  }
  if (endZ   != PETSC_NULL) {
    PetscValidPointer(endZ);
    *endZ   = mesh->locEndZ;
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshUpdateBoundingBox"
/*@
  MeshUpdateBoundingBox - This function updates the bounding box for the mesh based upon the current geometry.

  Not collective

  Input Parameter:
. mesh - The mesh

  Level: intermediate

.keywords: mesh, bounding box, update
.seealso: MeshSetBoundingBox(), MeshGetBoundingBox()
@*/
int MeshUpdateBoundingBox(Mesh mesh)
{
  MPI_Comm   comm;
  PetscTruth isPeriodic, isSetup;
  PetscReal  locStartX, locStartY, locStartZ;
  PetscReal  locEndX,   locEndY,   locEndZ;
  PetscReal  startX, startY, startZ;
  PetscReal  endX,   endY,   endZ;
  int        ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = PetscObjectGetComm((PetscObject) mesh, &comm);                                                   CHKERRQ(ierr);
  ierr = MeshIsPeriodic(mesh, &isPeriodic);                                                               CHKERRQ(ierr);
  if (isPeriodic == PETSC_TRUE) PetscFunctionReturn(0);
  /* Calculate the local bounding box */
  ierr = (*mesh->ops->updateboundingbox)(mesh);                                                           CHKERRQ(ierr);
  /* Calculate global bounding box */
  ierr = MeshGetLocalBoundingBox(mesh, &locStartX, &locStartY, &locStartZ, &locEndX, &locEndY, &locEndZ); CHKERRQ(ierr);
  ierr = MPI_Allreduce(&locStartX, &startX, 1, MPIU_REAL, MPI_MIN, comm);                                 CHKERRQ(ierr);
  ierr = MPI_Allreduce(&locEndX,   &endX,   1, MPIU_REAL, MPI_MAX, comm);                                 CHKERRQ(ierr);
  ierr = MPI_Allreduce(&locStartY, &startY, 1, MPIU_REAL, MPI_MIN, comm);                                 CHKERRQ(ierr);
  ierr = MPI_Allreduce(&locEndY,   &endY,   1, MPIU_REAL, MPI_MAX, comm);                                 CHKERRQ(ierr);
  ierr = MPI_Allreduce(&locStartZ, &startZ, 1, MPIU_REAL, MPI_MIN, comm);                                 CHKERRQ(ierr);
  ierr = MPI_Allreduce(&locEndZ,   &endZ,   1, MPIU_REAL, MPI_MAX, comm);                                 CHKERRQ(ierr);
  /* Pretend that the mesh is not setup */
  isSetup           = mesh->setupcalled;
  mesh->setupcalled = PETSC_FALSE;
  ierr = MeshSetBoundingBox(mesh, startX, startY, startZ, endX, endY, endZ);                              CHKERRQ(ierr);
  mesh->setupcalled = isSetup;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetPartition"
/*@
  MeshGetPartition - Returns the Partition object for this Mesh.

  Not collective

  Input Parameter:
. mesh - The mesh

  Output Parameter:
. part - The partition

  Level: intermediate

.keywords: Mesh, Partition, get
.seealso: PartitionGetMesh()
@*/
int MeshGetPartition(Mesh mesh, Partition *part)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(part);
  *part = mesh->part;
  PetscFunctionReturn(0);
}
#undef  __FUNCT__
#define __FUNCT__ "MeshSetMovement"
/*@
  MeshSetMovement - This function sets the dimension of the Mesh. This may only be done before a mesh is generated.

  Not collective

  Input Parameters:
+ mesh - The mesh
- dim  - The mesh dimension

  Level: intermediate

.keywords: Mesh, dimension
.seealso: MeshGetMovement(), MeshGetMover()
@*/
int MeshSetMovement(Mesh mesh, PetscTruth isMoving)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  mesh->isMoving = isMoving;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetMovement"
/*@
  MeshGetMovement - This function determines whether the Mesh is moving.

  Not collective

  Input Parameter:
. mesh     - The mesh

  Output Parameter:
. isMoving - The flag for mesh movement

  Level: intermediate

.keywords: Mesh, movement
.seealso: MeshSetMovement(), MeshGetMover()
@*/
int MeshGetMovement(Mesh mesh, PetscTruth *isMoving)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(isMoving);
  *isMoving = mesh->isMoving;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetMaxDegree"
/*@
  MeshGetMaxDegree - This function returns the maximum degree of any node in the mesh.

  Not collective

  Input Parameter:
. mesh      - The mesh

  Output Parameter:
. maxDegree - The maximum degree of any node

  Level: intermediate

.keywords: Mesh, degree
.seealso: MeshSetFromOptions()
@*/
int MeshGetMaxDegree(Mesh mesh, int *maxDegree)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(maxDegree);
  *maxDegree = mesh->maxDegree;
  PetscFunctionReturn(0);
}

/*------------------------------------------- Mesh Boundary Query Functions ------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "MeshGetNumBoundaries"
/*@
  MeshGetNumBoundaries - Returns the number of boundaries in the mesh.

  Not collective

  Input Parameter:
. mesh  - The mesh

  Output Parameter:
. numBd - The number boundaries

  Level: intermediate

.keywords: mesh, boundary
.seealso: MeshGetBoundaryStart(), MeshGetBoundarySize()
@*/
int MeshGetNumBoundaries(Mesh mesh, int *numBd) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(numBd);
  *numBd = mesh->numBd;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetBoundarySize"
/*@
  MeshGetBoundarySize - Returns the number of nodes on the specified boundary.

  Not collective

  Input Parameters:
+ mesh     - The mesh
- boundary - The boundary marker

  Output Parameter:
. size     - The number of nodes on the boundary

  Level: intermediate

.keywords: mesh, boundary
.seealso: MeshGetBoundaryStart()
@*/
int MeshGetBoundarySize(Mesh mesh, int boundary, int *size)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(size);
  ierr = (*mesh->ops->getboundarysize)(mesh, boundary, size);                                            CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetBoundaryIndex"
/*@
  MeshGetBoundaryIndex - Returns the index of the specified boundary.

  Not collective

  Input Parameters:
+ mesh     - The mesh
- boundary - The boundary marker

  Output Parameter:
. index - The index of the boundary

  Level: intermediate

.keywords: mesh, boundary
.seealso: MeshGetBoundarySize()
@*/
int MeshGetBoundaryIndex(Mesh mesh, int boundary, int *index)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(index);
  ierr = (*mesh->ops->getboundaryindex)(mesh, boundary, index);                                          CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetBoundaryStart"
/*@
  MeshGetBoundaryStart - Retrieves the canonical node number of the first node
  on the specified boundary.

  Not collective

  Input Parameters:
+ mesh     - The mesh
. boundary - The boundary marker
- ghost    - Flag for including ghost nodes

  Output Parameter:
. node     - The canonical node number

  Note:
$ This is typically used as an iteration construct with MeshGetBoundaryNext(),
$ for example:
$
$ MeshGetBoundaryStart(mesh, OUTER_BD, &node, PETSC_FALSE);
$ while (node >= 0) {
$   PetscPrintf(PETSC_COMM_SELF, "I have boundary node %d\n", node);
$   MeshGetBoundaryNext(mesh, OUTER_BD, &node, PETSC_FALSE);
$ }

  Level: intermediate

.keywords: mesh, boundary, iterator
.seealso: MeshGetBoundaryNext()
@*/
int MeshGetBoundaryStart(Mesh mesh, int boundary, PetscTruth ghost, int *node)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(node);
  ierr = (*mesh->ops->getboundarystart)(mesh, boundary, ghost, node);                                    CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetBoundaryNext"
/*@
  MeshGetBoundaryNext - Retrieves the canonical node number of the next node
  on the specified boundary, with -1 indicating boundary termination.

  Not collective

  Input Parameters:
+ mesh     - The mesh
. boundary - The boundary marker
- ghost    - Flag for including ghost nodes

  Output Parameter:
. node     - The canonical node number

  Note:
$ This is typically used as an iteration construct with MeshGetBoundaryStart(),
$ for example:
$
$ MeshGetBoundaryStart(mesh, OUTER_BD, &node, PETSC_FALSE);
$ while (node >= 0) {
$   PetscPrintf(PETSC_COMM_SELF, "I have boundary node %d\n", node);
$   MeshGetBoundaryNext(mesh, OUTER_BD, &node, PETSC_FALSE);
$ }

  Also, this only returns nodes which lie in the given domain, thus the above
  loop would run in parallel on all processors, returning different nodes on each.

  Level: intermediate

.keywords: mesh, boundary, iterator
.seealso: MeshGetBoundaryStart()
@*/
int MeshGetBoundaryNext(Mesh mesh, int boundary, PetscTruth ghost, int *node)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(node);
  ierr = (*mesh->ops->getboundarynext)(mesh, boundary, ghost, node);                                     CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetActiveBoundary"
/*@
  MeshGetActiveBoundary - Retrieves the boundary marker for the boundary
  last iterated over, or -1 if no iteration has taken place.

  Not collective

  Input Parameter:
. mesh     - The mesh

  Output Parameter:
. boundary - The boundary marker

  Level: advanced

.keywords: grid, boundary, iterator
.seealso: GridGetBoundaryNext(), MeshGetBoundaryStart()
@*/
int MeshGetActiveBoundary(Mesh mesh, int *boundary)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(boundary);
  ierr = (*mesh->ops->getactiveboundary)(mesh, boundary);                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*--------------------------------------------- Mesh Node Query Functions --------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "MeshGetNodeBoundary"
/*@
  MeshGetNodeBoundary - Returns the boundary on which the node lies, or 0 for interior nodes.

  Not collective

  Input Parameters:
+ mesh  - The mesh
- node  - The node

  Output Parameter:
. bd    - The boundary

  Level: intermediate

.keywords: mesh, boundary, node
.seealso: MeshGetBoundarySize()
@*/
int MeshGetNodeBoundary(Mesh mesh, int node, int *bd)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(bd);
  ierr = (*mesh->ops->getnodeboundary)(mesh, node, bd);                                                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshNodeIsVertex"
/*@
  MeshNodeIsVertex - This function determines whether or not a node is a vertex.

  Not collective

  Input Parameters:
+ mesh     - The mesh
- node     - The node

  Output Parameter:
. isVertex - The flag for a vertex

  Note:
  Vertices are nodes which are connected to edges.

  Level: intermediate

.keywords: mesh, node, vertex
.seealso: MeshGetNodeCoords()
@*/
int MeshNodeIsVertex(Mesh mesh, int node, PetscTruth *isVertex)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(isVertex);
  ierr = (*mesh->ops->nodeisvertex)(mesh, node, isVertex);                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNodeCoords"
/*@
  MeshGetNodeCoords - Retrieves the coordinates of a selected node

  Input Parameters:
+ mesh  - The mesh
- node  - The local node number

  Output Parameters:
. x,y,z - The coordinates

  Level: intermediate

.keywords: mesh, hole, coordinates
.seealso: MeshGetBoundaryStart()
@*/
int MeshGetNodeCoords(Mesh mesh, int node, double *x, double *y, double *z)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (x != PETSC_NULL) PetscValidScalarPointer(x);
  if (y != PETSC_NULL) PetscValidScalarPointer(y);
  if (z != PETSC_NULL) PetscValidScalarPointer(z);
  ierr = (*mesh->ops->getnodecoords)(mesh, node, x, y, z);                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSetNodeCoords"
/*@
  MeshSetNodeCoords - Sets the coordinates of a selected node

  Collective on Mesh

  Input Parameters:
+ mesh  - The mesh
. node  - The local node number
- x,y,z - The coordinates

  Level: intermediate

.keywords: mesh, node, coordinates
.seealso: MeshGetBoundaryStart()
@*/
int MeshSetNodeCoords(Mesh mesh, int node, double x, double y, double z)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = (*mesh->ops->setnodecoords)(mesh, node, x, y, z);                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNodeCoordsSaved"
/*@
  MeshGetNodeCoordsSaved - Retrieves the coordinates of a selected node
  from those saved with MeshSaveMesh().

  Input Parameters:
+ mesh  - The mesh
- node  - The canonical node number

  Output Parameters:
. x,y,z - The coordinates

  Level: intermediate

.keywords: mesh, node, coordinates
.seealso: MeshGetNodeCoords(), MeshSaveMesh()
@*/
int MeshGetNodeCoordsSaved(Mesh mesh, int node, double *x, double *y, double *z)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (x != PETSC_NULL) PetscValidScalarPointer(x);
  if (y != PETSC_NULL) PetscValidScalarPointer(y);
  if (z != PETSC_NULL) PetscValidScalarPointer(z);
  ierr = (*mesh->ops->getnodecoordssaved)(mesh, node, x, y, z);                                           CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNearestNode"
/*@
  MeshGetNearestNode - This function returns the node nearest to
  the given point, or -1 if the closest node is not contained in
  the local mesh.

  Not collective

  Input Parameters:
+ mesh    - The mesh
. x,y,z   - The node coordinates
- outside - A flag to allow points outside the domain

  Output Parameter:
. node    - The nearest node

  Note:
  The outside flag allows points outside the domain to be tested. If this flag
  is PETSC_FALSE, and (x,y,z) does not lie in the global domain, then an error
  will result.

  Level: beginner

.keywords: mesh, node, point location
.seealso MeshLocatePoint()
@*/
int MeshGetNearestNode(Mesh mesh, double x, double y, double z, PetscTruth outside, int *node)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(node);
  ierr = (*mesh->ops->nearestnode)(mesh, x, y, z, outside, node);                                         CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNearestBdNode"
/*@
  MeshGetNearestBdNode - This function returns the boundary node nearest to
  the given point, or -1 if the closest node is not contained in the local mesh.

  Not collective

  Input Parameters:
+ mesh  - The mesh
- x,y,z - The node coordinates

  Output Parameter:
. node  - The nearest boundary node

  Level: beginner

.keywords: mesh, node, point location
.seealso MeshGetNearestNode(), MeshLocatePoint()
@*/
int MeshGetNearestBdNode(Mesh mesh, double x, double y, double z, int *node)
{
  Partition p;
  double    minDist = PETSC_MAX;
  double    dist, nx, ny, nz;
  int       minNode, globalMinNode, allMinNode;
  int       numCorners, startNode, endNode;
  int       elem, corner, nearNode, bd;
  int       ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(node);
  ierr = MeshGetNumCorners(mesh, &numCorners);                                                            CHKERRQ(ierr);
  ierr = MeshGetPartition(mesh, &p);                                                                      CHKERRQ(ierr);
  ierr = PartitionGetStartNode(p, &startNode);                                                            CHKERRQ(ierr);
  ierr = PartitionGetEndNode(p, &endNode);                                                                CHKERRQ(ierr);
  ierr = MeshLocatePoint(mesh, x, y, z, &elem);                                                           CHKERRQ(ierr);
  if (elem >= 0) {
    /* Find first boundary node */
    for(corner = 0; corner < numCorners; corner++) {
      ierr = MeshGetNodeFromElement(mesh, elem, corner, &minNode);                                        CHKERRQ(ierr);
      ierr = MeshGetNodeBoundary(mesh, minNode, &bd);                                                     CHKERRQ(ierr);
      if (bd != 0) {
        ierr = MeshGetNodeCoords(mesh, minNode, &nx, &ny, &nz);                                           CHKERRQ(ierr);
        minDist = PetscSqr(MeshPeriodicDiffX(mesh, nx - x)) + PetscSqr(MeshPeriodicDiffY(mesh, ny - y));
        break;
      }
    }
    if (corner == numCorners) SETERRQ1(PETSC_ERR_ARG_WRONG, "Element %d has no node on a boundary", elem);
    /* Find closest node with field defined */
    for(corner = 1; corner < numCorners; corner++) {
      ierr = MeshGetNodeFromElement(mesh, elem, corner, &nearNode);                                       CHKERRQ(ierr);
      ierr = MeshGetNodeCoords(mesh, nearNode, &nx, &ny, &nz);                                            CHKERRQ(ierr);
      ierr = MeshGetNodeBoundary(mesh, nearNode, &bd);                                                    CHKERRQ(ierr);
      dist = PetscSqr(nx - x) + PetscSqr(ny - y);
      if ((bd != 0) && (dist < minDist)) {
        minDist = dist;
        minNode = nearNode;
      }
    }
    ierr = PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode);                                   CHKERRQ(ierr);
  } else {
    minNode       = -1;
    globalMinNode = -1;
  }

  /* The minimum node might still be a ghost node */
  ierr = MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, mesh->comm);                     CHKERRQ(ierr);
  if ((allMinNode >= startNode) && (allMinNode < endNode)) {
    *node = allMinNode - startNode;
  } else {
    *node = -1;
  }
  if (allMinNode < 0) PetscFunctionReturn(1);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNodeSupport"
/*@
  MeshGetNodeSupport - This function get the canonical element numbers of all elements
  within the support of a given basis function. A call to MeshRestoreNodeSupport() must
  be made before another call to this function.

  Not collective

  Input Parameters:
+ mesh    - The mesh
. node    - The node containing the basis function
- elem    - An element containing the node, -1 for a search

  Output Parameters:
+ degree  - The degree of the node
- support - A list of the elements in the support

  Note:
  This function currently only returns the elements containing
  any given node, so some basis functions will have a wider
  support than this definition.

  Level: intermediate

.keywords: mesh, support
.seealso: MeshRestoreNodeSupport()
@*/
int MeshGetNodeSupport(Mesh mesh, int node, int elem, int *degree, int **support)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(degree);
  PetscValidPointer(support);
  if (mesh->supportTaken == PETSC_TRUE) {
    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Missing call to MeshRestoreNodeSupport()");
  }
  mesh->supportTaken = PETSC_TRUE;
  ierr = (*mesh->ops->getnodesupport)(mesh, node, elem, degree, support);                                 CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshRestoreNodeSupport"
/*@
  MeshRestoreNodeSupport - This function returns the support array from MeshGetNodeSupport().

  Not collective

  Input Parameters:
+ mesh    - The mesh
. node    - The node containing the basis function
. elem    - An element containing the node, -1 for a search
. degree  - The degree pf the node
- support - A list of the elements in the support

  Level: intermediate

.keywords: mesh, support
.seealso: MeshGetNodeSupport()
@*/
int MeshRestoreNodeSupport(Mesh mesh, int node, int elem, int *degree, int **support)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (*support != mesh->support) {
    SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid support argument");
  }
  mesh->supportTaken = PETSC_FALSE;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNodeOrdering"
/*@
  MeshGetNodeOrdering - This function gets the AO which was used to reorder the mesh nodes
  before partitioning. This is sometimes used to bandwdith reduce the mesh graph.

  Not collective

  Input Parameter:
. mesh  - The mesh

  Output Parameter:
. order - The node ordering, or PETSC_NULL if none exists

  Level: intermediate

.keywords: Mesh, node, ordering, AO
.seealso: PartitionGetNodeOrdering()
@*/
int MeshGetNodeOrdering(Mesh mesh, AO *order)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(order);
  *order = mesh->nodeOrdering;
  PetscFunctionReturn(0);
}

/*------------------------------------------- Mesh Element Query Functions -------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "MeshGetElementNeighbor"
/*@
  MeshGetElementNeighbor - This function retrieves the local element number of the element, or neighbor, opposing a
  given corner.

  Not collective

  Input Parameters:
+ mesh     - The mesh
. elem     - The element
- corner   - The corner, or element node number

  Output Parameter:
. neighbor - The neighbor opposing corner

  Level: intermediate

.keywords: mesh, element, neighbor
.seealso: MeshGetNodeCoords()
@*/
int MeshGetElementNeighbor(Mesh mesh, int elem, int corner, int *neighbor)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(neighbor);
  ierr = (*mesh->ops->getelemneighbor)(mesh, elem, corner, neighbor);                                     CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshLocatePoint"
/*@
  MeshLocatePoint - This function returns the element containing
  the given point, or -1 if it is not contained in the local mesh.

  Not collective

  Input Parameters:
+ mesh           - The mesh
- x,y,z          - The node coordinates

  Output Parameter:
. containingElem - An element containing the node, -1 for failure

  Level: beginner

.keywords: mesh, point location
.seealso MeshGetNearestNode(), MeshGetNodeSupport()
@*/
int MeshLocatePoint(Mesh mesh, double x, double y, double z, int *containingElem)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidIntPointer(containingElem);
  ierr = PetscLogEventBegin(MESH_LocatePoint, mesh, 0, 0, 0);                                             CHKERRQ(ierr);
  ierr = (*mesh->ops->locatepoint)(mesh, x, y, z, containingElem);                                       CHKERRQ(ierr);
  ierr = PetscLogEventEnd(MESH_LocatePoint, mesh, 0, 0, 0);                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*--------------------------------------------- Mesh Hole Query Functions --------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "MeshSetHoleCoords"
/*@
  MeshSetHoleCoords - Sets the coordinates of holes in the mesh

  Collective on Mesh

  Input Parameter:
+ mesh   - The mesh
. num    - The number of holes
- coords - The coordinates

  Level: advanced

.keywords: mesh, hole, coordinates
.seealso: MeshGetBoundaryStart()
@*/
int MeshSetHoleCoords(Mesh mesh, int num, Vec coords)
{
  int          dim = mesh->dim;
  PetscScalar *array;
  int          size, vecDim;
  int          h, d;
  int          ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh,   MESH_COOKIE);
  PetscValidHeaderSpecific(coords, VEC_COOKIE);
  if (num != mesh->numHoles) {
    mesh->numHoles = num;
    if (mesh->holes != PETSC_NULL) {
      ierr = PetscFree(mesh->holes);                                                                      CHKERRQ(ierr);
    }
    if (mesh->numHoles > 0) {
      ierr = PetscMalloc(mesh->numHoles*dim * sizeof(double), &mesh->holes);                              CHKERRQ(ierr);
    } else {
      mesh->holes = PETSC_NULL;
    }
  }
  ierr = VecGetArray(coords, &array);                                                                     CHKERRQ(ierr);
  ierr = VecGetLocalSize(coords, &size);                                                                  CHKERRQ(ierr);
  if (size == mesh->numHoles*dim) {
    ierr = PetscMemcpy(mesh->holes, array, mesh->numHoles*dim * sizeof(double));                          CHKERRQ(ierr);
  } else if ((size > mesh->numHoles*dim) && (size%mesh->numHoles == 0)) {
    vecDim = (int) (size/mesh->numHoles);
    for(h = 0; h < mesh->numHoles; h++) {
      for(d = 0; d < dim; d++) {
        mesh->holes[h*dim+d] = PetscRealPart(array[h*vecDim+d]);
      }
    }
  } else {
    SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid vector of hole coordinates");
  }
  ierr = VecRestoreArray(coords, &array);                                                                 CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*------------------------------------------ Mesh Embedding Query Functions ------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "MeshGetElementFromNode"
/*@
  MeshGetElementFromNode - This function retrieves a local element number from the local
  node number. Notice that this could be nondeterministic.

  Not collective

  Input Parameters:
+ mesh - The mesh
- node - The number

  Output Parameter:
. elem - The element

  Level: intermediate

.keywords: mesh, node, element
.seealso: MeshGetNodeFromElement()
@*/
int MeshGetElementFromNode(Mesh mesh, int node, int *elem)
{
  Partition part;
  int       numElements, numCorners;
  int       e, c, n;
  int       ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(elem);
  ierr = MeshGetPartition(mesh, &part);                                                                  CHKERRQ(ierr);
  ierr = PartitionGetNumElements(part, &numElements);                                                    CHKERRQ(ierr);
  ierr = MeshGetNumCorners(mesh, &numCorners);                                                           CHKERRQ(ierr);
  for(e = 0; e < numElements; e++) {
    for(c = 0; c < numCorners; c++){
      ierr = MeshGetNodeFromElement(mesh, e, c, &n);                                                     CHKERRQ(ierr);
      if (n == node) {
        *elem = e;
        PetscFunctionReturn(0);
      }
    }
  }
  PetscFunctionReturn(-1);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetBdElementFromEdge"
/*@
  MeshGetBdElementFromEdge - This function retrieves the local element number of a
  boundary element from the local edge number.

  Not collective

  Input Parameters:
+ mesh - The mesh
- edge - The edge

  Output Parameter:
. elem - The element along the boundary

  Level: intermediate

.keywords: mesh, element, edge
.seealso: MeshGetNodeFromElement(), MeshGetNodeCoords()
@*/
int MeshGetBdElementFromEdge(Mesh mesh, int edge, int *elem)
{
#if 0
  int  neighbor; 
#endif
  int  degree;
  int *support;
  int  startNode,   endNode;
  int  startCorner, endCorner, numCorners;
  int  sElem, corner, node;
  int  ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  /* This should be a default function, which is only used if the implementation
     does not have an optimized implementation
  */
  /* Search for element containing edge */
  ierr = MeshGetNumCorners(mesh, &numCorners);                                                            CHKERRQ(ierr);
  ierr = MeshGetNodeFromEdge(mesh, edge, 0, &startNode);                                                  CHKERRQ(ierr);
  ierr = MeshGetNodeFromEdge(mesh, edge, 1, &endNode);                                                    CHKERRQ(ierr);
  /* Loop over nodes on each element in the support of the node */
  ierr = MeshGetNodeSupport(mesh, startNode, 0, &degree, &support);                                       CHKERRQ(ierr);
  for(sElem = 0; sElem < degree; sElem++) {
    for(corner = 0, startCorner = -1, endCorner = -1; corner < numCorners; corner++) {
      ierr = MeshGetNodeFromElement(mesh, support[sElem], corner, &node);                                 CHKERRQ(ierr);
      if (node == startNode) {
        startCorner = corner;
      } else if (node == endNode) {
        endCorner   = corner;
      }
    }
    if ((startCorner >= 0) && (endCorner >= 0)) break;
  }
  if (sElem == degree) SETERRQ1(PETSC_ERR_PLIB, "Edge %d not found in element list", edge);
  *elem = support[sElem];
  ierr = MeshRestoreNodeSupport(mesh, startNode, 0, &degree, &support);                                   CHKERRQ(ierr);
#if 0
  ierr = MeshGetNeighbor(mesh, *elem, ((startCorner+endCorner)*2)%3, &neighbor);                          CHKERRQ(ierr);  
  if (neighbor != -1)
    SETERRQ2(PETSC_ERR_ARG_WRONG, "Edge is not on a boundary since elem %d has neighbor %d", *elem, neighbor);
#endif
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNodeFromElement"
/*@
  MeshGetNodeFromElement - This function retrieves the local node number from the local
  element number and the corner.

  Not collective

  Input Parameters:
+ mesh   - The mesh
. elem   - The element
- corner - The corner, or element node number

  Output Parameter:
. node   - The local node number

  Level: intermediate

.keywords: mesh, node
.seealso: MeshGetNodeCoords()
@*/
int MeshGetNodeFromElement(Mesh mesh, int elem, int corner, int *node)
{
  int numCorners         = mesh->numCorners;
  int numOverlapElements = mesh->numFaces;
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(node);
  if (mesh->part != PETSC_NULL) {
    ierr  = PartitionGetNumOverlapElements(mesh->part, &numOverlapElements);                              CHKERRQ(ierr);
  }
  if ((elem < 0) || (elem >= numOverlapElements)) {
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node %d should be in [0,%d)", elem, numOverlapElements);
  }
  if ((corner < 0) || (corner >= numCorners)) {
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid corner %d should be in [0,%d)", elem, numCorners);
  }
  ierr = (*mesh->ops->getnodefromelement)(mesh, elem, corner, node);                                     CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetNodeFromEdge"
/*@
  MeshGetNodeFromEdge - This function retrieves the local node number from the local
  edge number and the corner.

  Not collective

  Input Parameters:
+ mesh   - The mesh
. edge   - The edge
- corner - The corner, or edge node number (0 or 1)

  Output Parameter:
. node   - The local node number

  Level: intermediate

.keywords: mesh, node, edge
.seealso: MeshGetNodeCoords()
@*/
int MeshGetNodeFromEdge(Mesh mesh, int edge, int corner, int *node)
{
  int numOverlapEdges = mesh->numEdges;
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(node);
  if ((corner < 0) || (corner >= 2)) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid corner %d should be in [0,2)", edge);
  if (mesh->part != PETSC_NULL) {
    ierr  = PartitionGetNumOverlapEdges(mesh->part, &numOverlapEdges);                                    CHKERRQ(ierr);
  }
  if ((edge < 0) || (edge >= numOverlapEdges)) {
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid edge %d should be in [0,%d)", edge, numOverlapEdges);
  }
  ierr = (*mesh->ops->getnodefromedge)(mesh, edge, corner, node);                                         CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetMidnodeFromEdge"
/*@
  MeshGetMidnodeFromEdge - This function retrieves the local node number of the midnode from the local edge number.

  Not collective

  Input Parameters:
+ mesh    - The mesh
- edge    - The edge

  Output Parameter:
. midnode - The local node number of the midnode

  Level: intermediate

.keywords: mesh, midnode, edge
.seealso: MeshGetNodeFromElement(), MeshGetNodeCoords()
@*/
int MeshGetMidnodeFromEdge(Mesh mesh, int edge, int *midnode)
{
  int  degree;
  int *support;
  int  startNode,   endNode;
  int  numCorners;
  int  startCorner = -1;
  int  endCorner   = -1;
  int  elem, corner, node;
  int  ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  /* This should be a default function, which is only used if the implementation
     does not have an optimized implementation
  */
  ierr = MeshGetNumCorners(mesh, &numCorners);                                                            CHKERRQ(ierr);
  if (numCorners == 3) PetscFunctionReturn(0);
  if (numCorners != 6) SETERRQ1(PETSC_ERR_ARG_WRONG, "Cannot find midnode for %d corners", numCorners);

  /* Search for midnode on edge */
  ierr = MeshGetNodeFromEdge(mesh, edge, 0, &startNode);                                                  CHKERRQ(ierr);
  ierr = MeshGetNodeFromEdge(mesh, edge, 1, &endNode);                                                    CHKERRQ(ierr);
  /* Loop over nodes on each element in the support of the node */
  ierr = MeshGetNodeSupport(mesh, startNode, 0, &degree, &support);                                       CHKERRQ(ierr);
  for(elem = 0; elem < degree; elem++) {
    for(corner = 0, startCorner = -1, endCorner = -1; corner < numCorners; corner++) {
      ierr = MeshGetNodeFromElement(mesh, support[elem], corner, &node);                                  CHKERRQ(ierr);
      if (node == startNode) {
        startCorner = corner;
      } else if (node == endNode) {
        endCorner   = corner;
      }
    }
    if ((startCorner >= 0) && (endCorner >= 0)) break;
  }
  if (elem == degree) {
    if (startCorner < 0) SETERRQ1(PETSC_ERR_PLIB, "Invalid support of node %d did not contain node", startNode);
    if (endCorner   < 0) {
      for(elem = 0; elem < mesh->numFaces; elem++) {
        for(corner = 0, startCorner = -1, endCorner = -1; corner < numCorners; corner++) {
          ierr = MeshGetNodeFromElement(mesh, elem, corner, &node);                                      CHKERRQ(ierr);
          if (node == startNode) {
            startCorner = corner;
          } else if (node == endNode) {
            endCorner   = corner;
          }
        }
        if ((startCorner >= 0) && (endCorner >= 0)) break;
      }
      if (elem == mesh->numFaces) {
        SETERRQ3(PETSC_ERR_PLIB, "Edge %d (%d,%d) not found in element list", edge, startNode, endNode);
      } else {
        SETERRQ5(PETSC_ERR_PLIB, "Edge %d (%d,%d) found in element %d not in support list of node %d",
                 edge, startNode, endNode, elem, startNode);
      }
    }
  }

  ierr = MeshGetNodeFromElement(mesh, support[elem], ((startCorner + endCorner)*2)%3 + 3, midnode);       CHKERRQ(ierr);
  ierr = MeshRestoreNodeSupport(mesh, startNode, 0, &degree, &support);                                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
}