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

/*
     Defines the interface to the mesh functions
*/

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

/* Logging support */
int MESH_COOKIE;
int MESH_Reform, MESH_IsDistorted, MESH_Partition, MESH_SetupBoundary, MESH_MoveMesh, MESH_CalcNodeVelocities;
int MESH_CalcNodeAccelerations, MESH_CreateLocalCSR, MESH_CreateFullCSR, MESH_CreateDualCSR, MESH_LocatePoint;

/*------------------------------------------------ Generic Operations ------------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "MeshSetUp"
/*@
  MeshSetUp - Set up any required internal data structures for a Mesh.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Level: beginner

.keywords: Mesh, setup
.seealso: MeshDestroy()
@*/
int MeshSetUp(Mesh mesh)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (mesh->setupcalled) PetscFunctionReturn(0);
  if (mesh->ops->setup) {
    ierr = (*mesh->ops->setup)(mesh);                                                                     CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MeshSetTypeFromOptions"
/*
  MeshSetTypeFromOptions - Sets the type of mesh generation from user options.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Level: intermediate

.keywords: Mesh, set, options, database, type
.seealso: MeshSetFromOptions(), MeshSetType()
*/
static int MeshSetTypeFromOptions(Mesh mesh)
{
  PetscTruth opt;
  char      *defaultType;
  char       typeName[256];
  int        dim;
  int        ierr;

  PetscFunctionBegin;
  ierr = MeshGetDimension(mesh, &dim);                                                                    CHKERRQ(ierr);
  if (mesh->dim == -1) {
    dim  = 2;
    if (mesh->dict != PETSC_NULL) {
      ierr = ParameterDictGetInteger(mesh->dict, "dim", &dim);                                            CHKERRQ(ierr);
    }
    ierr = MeshSetDimension(mesh, dim);                                                                   CHKERRQ(ierr);
  }
  if (mesh->type_name != PETSC_NULL) {
    defaultType = mesh->type_name;
  } else {
    switch(dim)
    {
    case 1:
      defaultType = MESH_TRIANGULAR_1D;
      break;
    case 2:
      defaultType = MESH_TRIANGULAR_2D;
      break;
    default:
      SETERRQ1(PETSC_ERR_SUP, "Mesh dimension %d is not supported", dim);
    }
  }

  if (!MeshRegisterAllCalled) {
    ierr = MeshRegisterAll(PETSC_NULL);                                                                   CHKERRQ(ierr);
  }
  ierr = PetscOptionsList("-mesh_type", "Mesh generation method"," MeshSetType", MeshList, defaultType, typeName, 256, &opt);
  CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = MeshSetType(mesh, typeName);                                                                   CHKERRQ(ierr);
  } else {
    ierr = MeshSetType(mesh, defaultType);                                                                CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MeshSetFromOptions"
/*@
  MeshSetFromOptions - Sets various Mesh parameters from user options.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Notes:  To see all options, run your program with the -help option, or consult the users manual.
          Must be called after MeshCreate() but before the Mesh is used.

  Level: intermediate

.keywords: Mesh, set, options, database
.seealso: MeshCreate(), MeshPrintHelp(), MeshSetOptionsPrefix()
@*/
int MeshSetFromOptions(Mesh mesh)
{
  Partition  part;
  PetscTruth opt;
  int        ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = PetscOptionsBegin(mesh->comm, mesh->prefix, "Mesh options", "Mesh");                             CHKERRQ(ierr); 

  /* Handle generic mesh options */
  ierr = PetscOptionsHasName(PETSC_NULL, "-help", &opt);                                                  CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = MeshPrintHelp(mesh);                                                                           CHKERRQ(ierr);
  }
  ierr = PetscOptionsGetInt(mesh->prefix, "-mesh_dim", &mesh->dim, &opt);                                 CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(mesh->prefix, "-mesh_max_aspect_ratio", &mesh->maxAspectRatio, &opt);        CHKERRQ(ierr);

  /* Handle mesh type options */
  ierr = MeshSetTypeFromOptions(mesh);                                                                    CHKERRQ(ierr);

  /* Handle specific mesh options */
  if (mesh->ops->setfromoptions != PETSC_NULL) {
    ierr = (*mesh->ops->setfromoptions)(mesh);                                                            CHKERRQ(ierr);
  }

  ierr = PetscOptionsEnd();                                                                               CHKERRQ(ierr);

  /* Handle subobject options */
  ierr = MeshGetPartition(mesh, &part);                                                                   CHKERRQ(ierr);
  ierr = PartitionSetFromOptions(part);                                                                   CHKERRQ(ierr);

  ierr = MeshViewFromOptions(mesh, mesh->name);                                                           CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshViewFromOptions"
/*@
  MeshViewFromOptions - This function visualizes the mesh based upon user options.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Level: intermediate

.keywords: Mesh, view, options, database
.seealso: MeshSetFromOptions(), MeshView()
@*/
int MeshViewFromOptions(Mesh mesh, char *title)
{
  PetscViewer viewer;
  PetscDraw   draw;
  PetscTruth  opt;
  char       *titleStr;
  char        typeName[1024];
  char        fileName[PETSC_MAX_PATH_LEN];
  int         len;
  int         ierr;

  PetscFunctionBegin;
  ierr = PetscOptionsHasName(mesh->prefix, "-mesh_view", &opt);                                           CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = PetscOptionsGetString(mesh->prefix, "-mesh_view", typeName, 1024, &opt);                       CHKERRQ(ierr);
    ierr = PetscStrlen(typeName, &len);                                                                   CHKERRQ(ierr);
    if (len > 0) {
      ierr = PetscViewerCreate(mesh->comm, &viewer);                                                      CHKERRQ(ierr);
      ierr = PetscViewerSetType(viewer, typeName);                                                        CHKERRQ(ierr);
      ierr = PetscOptionsGetString(mesh->prefix, "-mesh_view_file", fileName, 1024, &opt);                CHKERRQ(ierr);
      if (opt == PETSC_TRUE) {
        ierr = PetscViewerSetFilename(viewer, fileName);                                                  CHKERRQ(ierr);
      } else {
        ierr = PetscViewerSetFilename(viewer, mesh->name);                                                CHKERRQ(ierr);
      }
      ierr = MeshView(mesh, viewer);                                                                      CHKERRQ(ierr);
      ierr = PetscViewerFlush(viewer);                                                                    CHKERRQ(ierr);
      ierr = PetscViewerDestroy(viewer);                                                                  CHKERRQ(ierr);
    } else {
      ierr = MeshView(mesh, PETSC_NULL);                                                                  CHKERRQ(ierr);
    }
  }
  ierr = PetscOptionsHasName(mesh->prefix, "-mesh_view_draw", &opt);                                      CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = PetscViewerDrawOpen(mesh->comm, 0, 0, 0, 0, 300, 300, &viewer);                                CHKERRQ(ierr);
    ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);                                                      CHKERRQ(ierr);
    if (title != PETSC_NULL) {
      titleStr = title;
    } else {
      ierr = PetscObjectName((PetscObject) mesh);                                                         CHKERRQ(ierr) ;
      titleStr = mesh->name;
    }
    ierr = PetscDrawSetTitle(draw, titleStr);                                                             CHKERRQ(ierr);
    ierr = MeshView(mesh, viewer);                                                                        CHKERRQ(ierr);
    ierr = PetscViewerFlush(viewer);                                                                      CHKERRQ(ierr);
    ierr = PetscDrawPause(draw);                                                                          CHKERRQ(ierr);
    ierr = PetscViewerDestroy(viewer);                                                                    CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshView"
/*@
  MeshView - Views a mesh object.

  Collective on Mesh

  Input Parameters:
+ mesh   - The mesh
- viewer - The viewer with which to view the mesh

  Level: beginner

.keywords: mesh, view
.seealso: MeshDestroy(), PetscViewerDrawOpen()
@*/
int MeshView(Mesh mesh, PetscViewer viewer)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (viewer == PETSC_NULL) {
    viewer = PETSC_VIEWER_STDOUT_SELF;
  } else {
    PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE);
  }
  ierr = (*mesh->ops->view)(mesh, viewer);                                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPreCopy_Private"
/*
  MeshPreCopy_Private - Executes hook for generic pre-copying actions.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Level: developer

.keywords: Mesh, copy
.seealso: MeshPosCopy_Private(), MeshCopy(), MeshDuplicate(), MeshRefine(), MeshReform()
*/
int MeshPreCopy_Private(Mesh mesh) {
  int (*precopy)(Mesh);
  int   ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  /* Specific operations */
  ierr = PetscObjectQueryFunction((PetscObject) mesh, "PreCopy", (void (**)(void)) &precopy);             CHKERRQ(ierr);
  if (precopy != PETSC_NULL) {
    ierr = (*precopy)(mesh);                                                                              CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPostCopy_Private"
/*
  MeshPostCopy_Private - Executes hook for generic post-copying actions.

  Collective on Mesh

  Input Parameters:
+ mesh    - The mesh
- newMesh - The new mesh

  Level: developer

.keywords: Mesh, copy
.seealso: MeshPosCopy_Private(), MeshCopy(), MeshDuplicate(), MeshRefine(), MeshReform()
*/
int MeshPostCopy_Private(Mesh mesh, Mesh newMesh) {
  int      (*postcopy)(Mesh, Mesh);
  void      *ctx;
  PetscTruth isMoving;
  int        highlightElement;
  int        ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh,    MESH_COOKIE);
  PetscValidHeaderSpecific(newMesh, MESH_COOKIE);
  /* Generic operations */
  ierr = MeshGetUserContext(mesh, &ctx);                                                                  CHKERRQ(ierr);
  ierr = MeshSetUserContext(newMesh, ctx);                                                                CHKERRQ(ierr);
  ierr = MeshGetHighlightElement(mesh, &highlightElement);                                                CHKERRQ(ierr);
  ierr = MeshSetHighlightElement(newMesh, highlightElement);                                              CHKERRQ(ierr);
  ierr = MeshGetMovement(mesh, &isMoving);                                                                CHKERRQ(ierr);
  ierr = MeshSetMovement(newMesh, isMoving);                                                              CHKERRQ(ierr);
  /* Specific operations */
  ierr = PetscObjectQueryFunction((PetscObject) mesh, "PostCopy", (void (**)(void)) &postcopy);           CHKERRQ(ierr);
  if (postcopy != PETSC_NULL) {
    ierr = (*postcopy)(mesh, newMesh);                                                                    CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshCopy"
/*@
  MeshCopy - Copies all the data from one mesh to another.

  Collective on Mesh

  Input Parameter:
. mesh    - The mesh

  Output Parameter:
. newmesh - The identical new mesh

  Level: beginner

.keywords: Mesh, copy
.seealso: MeshDuplicate()
@*/
int MeshCopy(Mesh mesh, Mesh newmesh)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh,    MESH_COOKIE);
  PetscValidHeaderSpecific(newmesh, MESH_COOKIE);
  PetscFunctionBegin;
  if (mesh->ops->copy == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, " ");
  ierr = MeshPreCopy_Private(mesh);                                                                       CHKERRQ(ierr);
  ierr = (*mesh->ops->copy)(mesh, newmesh);                                                               CHKERRQ(ierr);
  ierr = MeshPostCopy_Private(mesh, newmesh);                                                             CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshDuplicate"
/*@
  MeshDuplicate - Duplicates the Mesh.

  Collective on Mesh

  Input Parameter:
. mesh    - The mesh

  Output Parameter:
. newmesh - The duplicate mesh

  Level: beginner

.keywords: mesh, duplicate, copy
.seealso: MeshCopy()
@*/
int MeshDuplicate(Mesh mesh, Mesh *newmesh)
{
  int   ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(newmesh);
  if (mesh->ops->duplicate == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, " ");
  ierr = MeshPreCopy_Private(mesh);                                                                       CHKERRQ(ierr);
  ierr = (*mesh->ops->duplicate)(mesh, newmesh);                                                          CHKERRQ(ierr);
  ierr = MeshPostCopy_Private(mesh, *newmesh);                                                            CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshDestroy"
/*@
  MeshDestroy - Destroys a mesh object.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Level: beginner

.keywords: mesh, destroy
.seealso MeshCreate()
@*/
int MeshDestroy(Mesh mesh)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (--mesh->refct > 0) PetscFunctionReturn(0);
  ierr = (*mesh->ops->destroy)(mesh);                                                                     CHKERRQ(ierr);
  if (mesh->holes != PETSC_NULL) {
    ierr = PetscFree(mesh->holes);                                                                        CHKERRQ(ierr);
  }
  PetscLogObjectDestroy(mesh);
  PetscHeaderDestroy(mesh);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSetUserContext"
/*@
  MeshSetUserContext - Sets the optional user-defined context for a mesh object.

  Collective on Mesh

  Input Parameters:
+ mesh   - The mesh
- usrCtx - The optional user context

  Level: intermediate

.keywords: mesh, set, context
.seealso: MeshGetUserContext(), GridSetMeshContext()
@*/
int MeshSetUserContext(Mesh mesh, void *usrCtx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  mesh->usr = usrCtx;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshGetUserContext"
/*@
  MeshGetUserContext - Gets the optional user-defined context for a mesh object.

  Not collective

  Input Parameters:
. mesh - The mesh

  Output Parameters:
. usrCtx   - The optional user context

  Level: intermediate

.keywords: mesh, set, context
.seealso: MeshSetUserContext(), GridGetMeshContext()
@*/
int MeshGetUserContext(Mesh mesh, void **usrCtx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(usrCtx);
  *usrCtx = mesh->usr;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MeshSetOptionsPrefix"
/*@C
  MeshSetOptionsPrefix - Sets the prefix used for searching for all Mesh options in the database.

  Not collective

  Input Parameters:
+ mesh   - The mesh
- prefix - The prefix to prepend to all option names

  Notes:
  A hyphen (-) must NOT be given at the beginning of the prefix name.
  The first character of all runtime options is AUTOMATICALLY the hyphen.

  Level: intermediate

.keywords: Mesh, set, options, prefix, database
.seealso: MeshGetOptionsPrefix(), MeshAppendOptionsPrefix(), MeshSetFromOptions()
@*/
int MeshSetOptionsPrefix(Mesh mesh, char *prefix)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh, prefix);                                         CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MeshAppendOptionsPrefix"
/*@C
  MeshAppendOptionsPrefix - Appends to the prefix used for searching for all Mesh options in the database.

  Not collective

  Input Parameters:
+ mesh   - The mesh
- prefix - The prefix to prepend to all option names

  Notes:
  A hyphen (-) must NOT be given at the beginning of the prefix name.
  The first character of all runtime options is AUTOMATICALLY the hyphen.

  Level: intermediate

.keywords: Mesh, append, options, prefix, database
.seealso: MeshGetOptionsPrefix(), MeshSetOptionsPrefix(), MeshSetFromOptions()
@*/
int MeshAppendOptionsPrefix(Mesh mesh, char *prefix)
{
  int ierr;
  
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = PetscObjectAppendOptionsPrefix((PetscObject) mesh, prefix);                                      CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MeshGetOptionsPrefix"
/*@
  MeshGetOptionsPrefix - Returns the prefix used for searching for all Mesh options in the database.

  Input Parameter:
. mesh   - The mesh

  Output Parameter:
. prefix - A pointer to the prefix string used

  Level: intermediate

.keywords: Mesh, get, options, prefix, database
.seealso: MeshSetOptionsPrefix(), MeshSetOptionsPrefix(), MeshSetFromOptions()
@*/
int MeshGetOptionsPrefix(Mesh mesh, char **prefix)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = PetscObjectGetOptionsPrefix((PetscObject) mesh, prefix);                                         CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "MeshPrintHelp"
/*@
  MeshPrintHelp - Prints all options for the Mesh.

  Input Parameter:
. mesh - The mesh

  Options Database Keys:
$  -help, -h

  Level: intermediate

.keywords: Mesh, help
.seealso: MeshSetFromOptions()
@*/
int MeshPrintHelp(Mesh mesh)
{
  char p[64];
  int  ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);

  ierr = PetscStrcpy(p, "-");                                                                             CHKERRQ(ierr);
  if (mesh->prefix != PETSC_NULL) {
    ierr = PetscStrcat(p, mesh->prefix);                                                                  CHKERRQ(ierr);
  }

  (*PetscHelpPrintf)(mesh->comm, "Mesh options ------------------------------------------------\n");
  (*PetscHelpPrintf)(mesh->comm,"   %smesh_type <typename> : Sets the mesh type\n", p);
  (*PetscHelpPrintf)(mesh->comm,"   %smesh_dim  <num>      : Sets the dimension of the mesh\n", p);
  PetscFunctionReturn(0);
}

/*--------------------------------------------- Mesh-Specific Operations ---------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "MeshCheckBoundary_Triangular_1D"
/*
  MeshCheckBoundary_Triangular_1D - Checks a boundary for validity.

  Collective on Mesh

  Input Parameter:
. mesh        - The mesh

  Input Parameters from bdCtx:
+ numBD       - The number of closed boundaries in the geometry, or different markers
. numVertices - The umber of boundary points
. vertices    - The (x,y) coordinates of the boundary points
. markers     - The boundary markers for nodes, 0 indicates an interior point, each boundary must have a different marker

  Level: advanced

.keywords: mesh, boundary
.seealso: MeshSetBoundary(), MeshSetReformBoundary()
*/
static int MeshCheckBoundary_Triangular_1D(Mesh mesh, MeshBoundary2D *bdCtx) {
  int rank;
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = MPI_Comm_rank(mesh->comm, &rank);                                                                CHKERRQ(ierr);
  if (rank == 0) {
    if (bdCtx->numBd != bdCtx->numVertices) {
      SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "All boundaries must be points: %d != %d", bdCtx->numBd, bdCtx->numVertices);
    }
    if ((bdCtx->numVertices != 0) && (bdCtx->numVertices != 2)) {
      SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of boundary points %d", bdCtx->numVertices);
    }
    if (bdCtx->numVertices > 0) {
      PetscValidScalarPointer(bdCtx->vertices);
      PetscValidIntPointer(bdCtx->markers);
    }
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshCheckBoundary_Triangular_2D"
/*
  MeshCheckBoundary_Triangular_2D - Checks a boundary for validity.

  Collective on Mesh

  Input Parameter:
. mesh        - The mesh

  Input Parameters from bdCtx:
+ numBD       - The number of closed boundaries in the geometry, or different markers
. numVertices - The umber of boundary points
. vertices    - The (x,y) coordinates of the boundary points
. markers     - The boundary markers for nodes, 0 indicates an interior point, each boundary must have a different marker
. numSegments - The number of boundary segments
. segments    - The endpoints of boundary segments or PETSC_NULL
. segMarkers  - The boundary markers for each segment
. numHoles    - The number of holes
- holes       - The (x,y) coordinates of holes or PETSC_NULL

  Level: advanced

.keywords: mesh, boundary
.seealso: MeshSetBoundary(), MeshSetReformBoundary()
*/
static int MeshCheckBoundary_Triangular_2D(Mesh mesh, MeshBoundary2D *bdCtx) {
  int rank;
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = MPI_Comm_rank(mesh->comm, &rank);                                                                CHKERRQ(ierr);
  if (rank == 0) {
    if (bdCtx->numVertices < 3) {
      SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Insufficient number of boundary points %d", bdCtx->numVertices);
    }
    if (bdCtx->numSegments < 3) {
      SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of boundary segments %d", bdCtx->numSegments);
    }
    if (bdCtx->numHoles    < 0) {
      SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of holes %d", bdCtx->numHoles);
    }
    if (bdCtx->numVertices > 0) {
      PetscValidScalarPointer(bdCtx->vertices);
      PetscValidIntPointer(bdCtx->markers);
    }
    if (bdCtx->numSegments > 0) {
      PetscValidIntPointer(bdCtx->segments);
      PetscValidIntPointer(bdCtx->segMarkers);
    }
    if (bdCtx->numHoles    > 0) {
      PetscValidScalarPointer(bdCtx->holes);
    }
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshCheckBoundary"
/*@
  MeshCheckBoundary - Checks the mesh boundary for validity.

  Collective on Mesh

  Input Parameters:
+ mesh  - The mesh
- bdCtx - The MeshBoundary

  Level: intermediate

.keywords: mesh, partition
.seealso: MeshSetBoundary(), MeshSetReformBoundary()
@*/
int MeshCheckBoundary(Mesh mesh, MeshBoundary2D *bdCtx) {
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  /* This needs to go in MeshBoundary later */
  if (mesh->dim == 1) {
    ierr = MeshCheckBoundary_Triangular_1D(mesh, bdCtx);                                                  CHKERRQ(ierr);
  } else if (mesh->dim == 2) {
    ierr = MeshCheckBoundary_Triangular_2D(mesh, bdCtx);                                                  CHKERRQ(ierr);
  } else {
    SETERRQ(PETSC_ERR_ARG_WRONG, "Mehs dimension has to be set in order to validate boundary");
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSetBoundary"
/*@
  MeshSetBoundary - Store a boundary to use for mesh generation.

  Not collective

  Input Parameters:
+ mesh  - The mesh
- bdCtx - The MeshBoundary

  Level: intermediate

.keywords: mesh, boundary, refinement
.seealso: MeshReform()
@*/
int MeshSetBoundary(Mesh mesh, MeshBoundary2D *bdCtx) {
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(bdCtx);
  ierr = MeshCheckBoundary(mesh, bdCtx);                                                                  CHKERRQ(ierr);
  mesh->bdCtx = bdCtx;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSetReformBoundary"
/*@
  MeshSetReformBoundary - Store an alternate boundary to use for reforming.

  Not collective

  Input Parameters:
+ mesh  - The mesh
- bdCtx - The MeshBoundary

  Level: intermediate

.keywords: mesh, boundary, refinement
.seealso: MeshReform()
@*/
int MeshSetReformBoundary(Mesh mesh, MeshBoundary2D *bdCtx) {
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(bdCtx);
  ierr = MeshCheckBoundary(mesh, bdCtx);                                                                  CHKERRQ(ierr);
  mesh->bdCtxNew = bdCtx;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshPartition"
/*@
  MeshPartition - Partitions the mesh among the processors

  Collective on Mesh

  Input Parameter:
. mesh - the mesh

  Level: beginner

.keywords: mesh, partition
.seealso: MeshGetBoundaryStart()
@*/
int MeshPartition(Mesh mesh)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  ierr = PetscLogEventBegin(MESH_Partition, mesh, 0, 0, 0);                                                CHKERRQ(ierr);
  ierr = (*mesh->ops->partition)(mesh);                                                                   CHKERRQ(ierr);
  ierr = PetscLogEventEnd(MESH_Partition, mesh, 0, 0, 0);                                                  CHKERRQ(ierr);
  ierr = MeshUpdateBoundingBox(mesh);                                                                     CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
 
#undef  __FUNCT__
#define __FUNCT__ "MeshCoarsen"
/*@
  MeshCoarsen - Coarsen a mesh based on area constraints.

  Collective on Mesh

  Input Parameters:
+ mesh    - The initial mesh
- area    - A function which gives an area constraint when evaluated inside an element

  Output Parameter:
. newmesh - The coarse mesh

  Note:
  If PETSC_NULL is used for the 'area' argument, then the mesh consisting only of vertices
  is returned.

  Level: beginner

.keywords: mesh, coarsening
.seealso: MeshRefine(), MeshDestroy()
@*/
int MeshCoarsen(Mesh mesh, PointFunction area, Mesh *newmesh)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(newmesh);
  if (mesh->ops->coarsen == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, " ");
  ierr = MeshPreCopy_Private(mesh);                                                                       CHKERRQ(ierr);
  ierr = (*mesh->ops->coarsen)(mesh, area, newmesh);                                                      CHKERRQ(ierr);
  ierr = MeshPostCopy_Private(mesh, *newmesh);                                                            CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
 
#undef  __FUNCT__
#define __FUNCT__ "MeshRefine"
/*@
  MeshRefine - Refine a mesh based on area constraints.

  Input Parameters:
+ mesh    - The initial mesh
- area    - A function which gives an area constraint when evaluated inside an element

  Output Parameter:
. newmesh - The refined mesh

  Level: beginner

.keywords: mesh, refinement
.seealso: MeshCoarsen(), MeshDestroy()
@*/
int MeshRefine(Mesh mesh, PointFunction area, Mesh *newmesh)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(newmesh);
  if (mesh->ops->refine == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, " ");
  ierr = MeshPreCopy_Private(mesh);                                                                       CHKERRQ(ierr);
  ierr = (*mesh->ops->refine)(mesh, area, newmesh);                                                       CHKERRQ(ierr);
  ierr = MeshPostCopy_Private(mesh, *newmesh);                                                            CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshResetNodes"
/*@
  MeshResetNodes - Using the vertex and edge structure, move auxilliary nodes back to
  their default positions.

  Input Parameters:
+ mesh    - The mesh
- resetBd - The flag which indicates whether boundaries should also be reset

  Note:
  This function was developed to allow midnodes to be moved back onto the edges
  with which they were associated. Midnodes on edges between nodes on the same
  boundary are only reset if resetBd == PETSC_TRUE (to allow for curved boundaries).

  Level: advanced

.keywords: mesh, reset, node
.seealso: MeshReform(), MeshSetBounday()
@*/
int MeshResetNodes(Mesh mesh, PetscTruth resetBd)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  if (!mesh->ops->resetnodes) {
    SETERRQ(PETSC_ERR_SUP, "Not supported by this Mesh object");
  }
  ierr = (*mesh->ops->resetnodes)(mesh, resetBd);                                                         CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "MeshSaveMesh"
/*@
  MeshSaveMesh - This function saves the mesh coordinates.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Level: advanced

  Note:
  This function is meant to be used in conjuction with MeshMoveMesh(),
  and MeshRestoreMesh() so that the mesh may be moved, checked for
  distortion, and if necessary moved back.

.keywords: mesh, movement
.seealso: MeshMoveMesh(), MeshRestoreMesh()
@*/
int MeshSaveMesh(Mesh mesh) {
  int ierr;

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

#undef  __FUNCT__
#define __FUNCT__ "MeshRestoreMesh"
/*@
  MeshRestoreMesh - This function restores the mesh coordinates which were
  previously saved.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Level: advanced

  Note:
  This function is meant to be used in conjuction with MeshMoveMesh(),
  and MeshRestoreMesh() so that the mesh may be moved, checked for
  distortion, and if necessary moved back.

.keywords: mesh, movement
.seealso: MeshMoveMesh(), MeshSaveMesh()
@*/
int MeshRestoreMesh(Mesh mesh) {
  int ierr;

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

#undef  __FUNCT__
#define __FUNCT__ "MeshIsDistorted"
/*@
  MeshIsDistorted - This function checks for more than the maximum
  level of distortion in the mesh. It also returns an error when
  encountering elements with negative area.

  Collective on Mesh

  Input Parameter:
. mesh - The mesh

  Output Parameter:
. flag - Signals a distorted mesh

  Level: intermediate

.keywords: mesh, distortion, movement
.seealso: MeshMoveMesh()
@*/
int MeshIsDistorted(Mesh mesh, PetscTruth *flag) {
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mesh, MESH_COOKIE);
  PetscValidPointer(flag);
  ierr = PetscLogEventBegin(MESH_IsDistorted, mesh, 0, 0, 0);                                              CHKERRQ(ierr);
  /* Do not use CHKERRQ since we want the return value at the top level */
  ierr = (*mesh->ops->isdistorted)(mesh, flag);
  ierr = PetscLogEventEnd(MESH_IsDistorted, mesh, 0, 0, 0);                                                CHKERRQ(ierr);
  PetscFunctionReturn(ierr);
}