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

#include "src/grid/gridimpl.h"    /*I "grid.h" I*/

/* Loggging support */
int VAR_ORDER_COOKIE;

/*------------------------------------------------- Variable Ordering ------------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "VarOrderingCreate"
/*@C
  VarOrderingCreate - This function creates a global variable ordering.

  Collective on Grid

  Input Parameter:
. grid  - The underlying grid for the ordering

  Output Parameter:
. order - The ordering

  Level: beginner

.keywords: variable ordering, create
.seealso: VarOrderingCreateConstrained(), VarOrderingDestroy()
@*/
int VarOrderingCreate(Grid grid, VarOrdering *order)
{ 
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidPointer(order);
  ierr = GridSetUp(grid);                                                                                CHKERRQ(ierr);
  ierr = (*grid->ops->gridcreatevarordering)(grid, grid->cm, order);                                     CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingCreateConstrained"
/*@C
  VarOrderingCreateConstrained - This function creates a global variable ordering on the constraint space.

  Collective on Grid

  Input Parameter:
. grid  - The underlying grid for the ordering

  Output Parameter:
. order - The ordering

  Level: beginner

.keywords: variable ordering, create, constraint
.seealso: VarOrderingCreate(), VarOrderingDestroy()
@*/
int VarOrderingCreateConstrained(Grid grid, VarOrdering *order)
{
  FieldClassMap cm;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidPointer(order);
  ierr = GridSetUp(grid);                                                                                CHKERRQ(ierr);
  ierr = GridGetClassMap(grid, &cm);                                                                     CHKERRQ(ierr);
  ierr = (*grid->ops->gridcreatevarordering)(grid, cm, order);                                           CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingCreateGeneral"
/*@C
  VarOrderingCreateGeneral - This function creates a global variable ordering on the specified fields.

  Collective on Grid

  Input Parameters:
+ grid      - The underlying grid for the ordering
. numFields - The number of fields
- fields    - The fields

  Output Parameter:
. order     - The ordering

  Level: beginner

.keywords: variable ordering, create, constraint
.seealso: VarOrderingCreate(), VarOrderingDestroy()
@*/
int VarOrderingCreateGeneral(Grid grid, int numFields, int *fields, VarOrdering *order)
{
  FieldClassMap map, tempMap;
  PetscTruth    isConstrained, reduceSystem;
  int           numTotalFields, f;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidPointer(order);
  ierr = GridGetNumFields(grid, &numTotalFields);                                                         CHKERRQ(ierr);
  if (numFields > numTotalFields) {
    SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid number  %d of fields, grid only has %d", numFields, numTotalFields);
  }
  for(f = 0; f < numFields; f++) {
    if ((fields[f] < 0) || (fields[f] >= numTotalFields)) {
      SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", fields[f]);
    }
  }
  ierr = GridSetUp(grid);                                                                                 CHKERRQ(ierr);
  ierr = FieldClassMapCreateTriangular2D(grid, numFields, fields, &tempMap);                              CHKERRQ(ierr);
  ierr = GridIsConstrained(grid, &isConstrained);                                                         CHKERRQ(ierr);
  ierr = GridGetReduceSystem(grid, &reduceSystem);                                                        CHKERRQ(ierr);
  if (reduceSystem == PETSC_TRUE) {
    ierr = FieldClassMapConstrain(tempMap, grid, reduceSystem, isConstrained, &map);                      CHKERRQ(ierr);
    ierr = FieldClassMapDestroy(tempMap);                                                                 CHKERRQ(ierr);
  } else {
    map  = tempMap;
  }
  ierr = (*grid->ops->gridcreatevarordering)(grid, map, order);                                           CHKERRQ(ierr);
  ierr = FieldClassMapDestroy(map);                                                                       CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingCreateSubset"
/*@C
  VarOrderingCreateSubset - This function creates a new global variable ordering from an existing
  one on a subset of the fields.

  Collective on Grid

  Input Parameters:
+ order     - The original ordering
. numFields - The number of fields in the new ordering
. fields    - The fields in the new ordering
- contract  - The flag for contracting the indices

  Output Parameter:
. newOrder  - The new ordering

  Note:
  If contract is PETSC_FALSE, this function returns a copy of the original
  ordering, but with some fields no longer active.

  Level: beginner

.keywords: variable ordering, create
.seealso: VarOrderingCreate(), VarOrderingDestroy()
@*/
int VarOrderingCreateSubset(VarOrdering order, int numFields, int *fields, PetscTruth contract, VarOrdering *newOrder)
{
  FieldClassMap map, cm;
  int           f, g;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  PetscValidPointer(newOrder);
  ierr = VarOrderingGetClassMap(order, &map);                                                            CHKERRQ(ierr);
  /* Check for consistency */
  for(f = 0; f < numFields; f++) {
    for(g = 0; g < map->numFields; g++) {
      if (map->fields[g] == fields[f])
        break;
    }
    if (g == map->numFields) {
      SETERRQ1(PETSC_ERR_ARG_INCOMP, "Fields not a subset: field %d not in the existing order", fields[f]);
    }
  }
  if (contract == PETSC_TRUE) {
    SETERRQ(PETSC_ERR_SUP, "Coming soon");
  } else {
    ierr = FieldClassMapCreateSubset(map, numFields, fields, &cm);                                       CHKERRQ(ierr);
    ierr = VarOrderingDuplicate(order, newOrder);                                                        CHKERRQ(ierr);
    ierr = PetscObjectCompose((PetscObject) *newOrder, "ClassMap", (PetscObject) cm);                    CHKERRQ(ierr);
    ierr = FieldClassMapDestroy(cm);                                                                     CHKERRQ(ierr);
    /* Free extra localStart entries */
    for(g = 0; g < map->numFields; g++) {
      for(f = 0; f < numFields; f++) {
        if (map->fields[g] == fields[f])
          break;
      }
      if (f == numFields) {
        ierr = PetscFree((*newOrder)->localStart[map->fields[g]]);                                       CHKERRQ(ierr);
      }
    }
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingConstrain"
/*@C
  VarOrderingConstrain - This function constrains a previous variable ordering using the constraints
  already defined in the grid.

  Collective on VarOrdering

  Input Parameters:
+ grid     - The underlying grid for the ordering
- oldOrder - The original ordering

  Output Parameter:
. order    - The constrained ordering

  Level: beginner

.keywords variable ordering, finite element
.seealso VarOrderingDestroy()
@*/
int VarOrderingConstrain(Grid grid, VarOrdering oldOrder, VarOrdering *order)
{
  PetscTruth    isConstrained;
  FieldClassMap cm;
  int           numFields;
  int           field;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidPointer(order);
  ierr = GridSetUp(grid);                                                                                 CHKERRQ(ierr);
  ierr = GridIsConstrained(grid, &isConstrained);                                                         CHKERRQ(ierr);
  if (isConstrained == PETSC_TRUE) {
    ierr = GridGetNumFields(grid, &numFields);                                                            CHKERRQ(ierr);
    ierr = GridGetClassMap(grid, &cm);                                                                    CHKERRQ(ierr);
    for(field = 0; field < numFields; field++) {
      if (grid->fields[field].isConstrained == PETSC_TRUE) break;
    }
    if (field == numFields) SETERRQ(PETSC_ERR_ARG_WRONG, "No constrained fields");
    ierr = (*grid->ops->gridvarorderingconstrain)(grid, field, grid->constraintCtx, cm, oldOrder, order); CHKERRQ(ierr);
  } else {
    ierr = VarOrderingDuplicate(oldOrder, order);                                                         CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingCreateReduce"
/*@C
  VarOrderingCreateReduce - This function creates a global variable ordering on the reduction space.

  Collective on Grid

  Input Parameter:
. grid  - The underlying grid for the ordering

  Output Parameter:
. order - The reduction ordering

  Level: beginner

.keywords: variable ordering, create, reduction
.seealso: VarOrderingCreate(), VarOrderingDestroy()
@*/
int VarOrderingCreateReduce(Grid grid, VarOrdering *order)
{ 
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidPointer(order);
  ierr = GridSetUp(grid);                                                                                CHKERRQ(ierr);
  if (grid->reduceSystem == PETSC_TRUE) {
    ierr = (*grid->ops->gridcreatevarordering)(grid, grid->reductionCM, order);                          CHKERRQ(ierr);
  } else {
    ierr = (*grid->ops->gridcreatevarordering)(grid, grid->cm,          order);                          CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingCreateReduceGeneral"
/*@C
  VarOrderingCreateReduceGeneral - This function creates a global variable ordering on the reduction space.

  Collective on Grid

  Input Parameters:
+ grid        - The underlying grid for the ordering
. cm          - The class map
- reductionCM - The class map for the reduction space

  Output Parameter:
. order       - The reduction ordering

  Level: beginner

.keywords: variable ordering, create, reduction
.seealso: VarOrderingCreate(), VarOrderingDestroy()
@*/
int VarOrderingCreateReduceGeneral(Grid grid, FieldClassMap cm, FieldClassMap reductionCM, VarOrdering *order)
{ 
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid,        GRID_COOKIE);
  PetscValidHeaderSpecific(cm,          CLASS_MAP_COOKIE);
  PetscValidHeaderSpecific(reductionCM, CLASS_MAP_COOKIE);
  PetscValidPointer(order);
  ierr = GridSetUp(grid);                                                                                 CHKERRQ(ierr);
  if (grid->reduceSystem == PETSC_TRUE) {
    ierr = (*grid->ops->gridcreatevarordering)(grid, reductionCM, order);                                 CHKERRQ(ierr);
  } else {
    ierr = (*grid->ops->gridcreatevarordering)(grid, cm,          order);                                 CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingDestroy"
/*@C
  VarOrderingDestroy - This function destroys a global variable ordering.

  Not collective

  Input Parameter:
. order - The ordering

  Level: beginner

.keywords: variable ordering
.seealso: VarOrderingCreate()
@*/
int VarOrderingDestroy(VarOrdering order)
{
  FieldClassMap map;
  int           f;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  if (--order->refct > 0)
    PetscFunctionReturn(0);
  ierr = PetscObjectQuery((PetscObject) order, "ClassMap", (PetscObject *) &map);                        CHKERRQ(ierr);
  ierr = PetscFree(order->firstVar);                                                                     CHKERRQ(ierr);
  ierr = PetscFree(order->offsets);                                                                      CHKERRQ(ierr);
  if (order->localOffsets) {
    ierr = PetscFree(order->localOffsets);                                                               CHKERRQ(ierr);
  }
  for(f = 0; f < map->numFields; f++) {
    ierr = PetscFree(order->localStart[map->fields[f]]);                                                 CHKERRQ(ierr);
  }
  ierr = PetscFree(order->localStart);                                                                   CHKERRQ(ierr);
  PetscLogObjectDestroy(order);
  PetscHeaderDestroy(order);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingGetClassMap"
/*@
  VarOrderingGetClassMap - This function returns the class map for the variable ordering.

  Not collective

  Input Parameter:
. order - The variable ordering

  Output Parameter:
. map   - The field class map

  Level: intermediate

.keywords: variable ordering, class map
.seealso: VarOrderingGetLocalSize()
@*/
int VarOrderingGetClassMap(VarOrdering order, FieldClassMap *map)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  PetscValidPointer(map);
  ierr = PetscObjectQuery((PetscObject) order, "ClassMap", (PetscObject *) map);                         CHKERRQ(ierr);
  PetscValidHeaderSpecific(*map,  CLASS_MAP_COOKIE);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingGetNumTotalFields"
/*@
  VarOrderingGetNumTotalFields - Gets the total number of fields in the associated grid.

  Not collective

  Input Parameter:
. order     - The variable ordering

  Output Parameter:
. numFields - The total number fields

  Level: intermediate

.keywords: variable ordering, grid, field
.seealso: VarOrderingGetSize(), VarOrderingGetLocalSize()
@*/
int VarOrderingGetNumTotalFields(VarOrdering order, int *numFields)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  PetscValidPointer(numFields);
  *numFields = order->numTotalFields;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingGetSize"
/*@
  VarOrderingGetSize - Gets the number of global variables in a variable ordering.

  Not collective

  Input Parameter:
. order - The variable ordering

  Output Parameter:
. size  - The number of global variables

  Level: intermediate

.keywords: variable ordering, grid
.seealso: VarOrderingGetLocalSize()
@*/
int VarOrderingGetSize(VarOrdering order, int *size)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  PetscValidPointer(size);
  *size = order->numVars;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingGetLocalSize"
/*@
  VarOrderingGetLocalSize - Gets the number of local variables in a variable ordering.

  Not collective

  Input Parameter:
. order - The variable ordering

  Output Parameter:
. size  - The number of local variables

  Level: intermediate

.keywords: variable ordering, grid
.seealso: VarOrderingGetSize()
@*/
int VarOrderingGetLocalSize(VarOrdering order, int *size)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  PetscValidPointer(size);
  *size = order->numLocVars;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingGetFirst"
/*@
  VarOrderingGetFirst - Gets the layout of variables across domains.

  Not collective

  Input Parameter:
. order - The variable ordering

  Output Parameter:
. first - An array of the first variable in each domain, and the
          total number of variables at the end

  Level: intermediate

.keywords: variable ordering, grid
.seealso: VarOrderingGetSize()
@*/
int VarOrderingGetFirst(VarOrdering order, int **first)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  PetscValidPointer(first);
  *first = order->firstVar;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingSerialize"
/*@ 
  VarOrderingSerialize - This function stores or recreates a variable ordering using a viewer for a binary file.

  Collective on Grid

  Input Parameter:
+ map    - The class map
. viewer - The viewer context
- store  - This flag is PETSC_TRUE is data is being written, otherwise it will be read

  Output Parameter:
. order  - The ordering

  Level: beginner

.keywords: variable ordering, serialize
.seealso: GridSerialize()
@*/
int VarOrderingSerialize(FieldClassMap map, VarOrdering *order, PetscViewer viewer, PetscTruth store)
{
  VarOrdering o;
  int         fd;
  int         cookie;
  int         numOverlapNodes = map->numOverlapNodes;
  int         numGhostNodes   = map->numGhostNodes;
  int         numClasses      = map->numClasses;
  int         zero            = 0;
  int         one             = 1;
  int         numProcs, hasLocOffsets;
  int         f, field;
  PetscTruth  match;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE);
  PetscValidPointer(order);

  ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &match);                             CHKERRQ(ierr);
  if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONG, "Must be binary viewer");
  ierr = PetscViewerBinaryGetDescriptor(viewer, &fd);                                                     CHKERRQ(ierr);
  ierr = MPI_Comm_size(map->comm, &numProcs);                                                             CHKERRQ(ierr);

  if (store) {
    o = *order;
    PetscValidHeaderSpecific(o, VAR_ORDER_COOKIE);
    ierr = PetscBinaryWrite(fd, &o->cookie,             1,               PETSC_INT, 0);                   CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &o->numVars,            1,               PETSC_INT, 0);                   CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &o->numLocVars,         1,               PETSC_INT, 0);                   CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &o->numOverlapVars,     1,               PETSC_INT, 0);                   CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &o->numNewVars,         1,               PETSC_INT, 0);                   CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &o->numLocNewVars,      1,               PETSC_INT, 0);                   CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &o->numTotalFields,     1,               PETSC_INT, 0);                   CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  o->firstVar,           numProcs+1,      PETSC_INT, 0);                   CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  o->offsets,            numOverlapNodes, PETSC_INT, 0);                   CHKERRQ(ierr);
    if (o->localOffsets == PETSC_NULL) {
      ierr = PetscBinaryWrite(fd, &zero,                1,               PETSC_INT, 0);                   CHKERRQ(ierr);
    } else {
      ierr = PetscBinaryWrite(fd, &one,                 1,               PETSC_INT, 0);                   CHKERRQ(ierr);
      ierr = PetscBinaryWrite(fd,  o->localOffsets,     numGhostNodes,   PETSC_INT, 0);                   CHKERRQ(ierr);
    }
    for(f = 0; f < map->numFields; f++) {
      field = map->fields[f];
      ierr = PetscBinaryWrite(fd, o->localStart[field], numClasses,      PETSC_INT, 0);                   CHKERRQ(ierr);
    }
  } else {
    ierr = PetscBinaryRead(fd, &cookie,                1,               PETSC_INT);                       CHKERRQ(ierr);
    if (cookie != VAR_ORDER_COOKIE) SETERRQ(PETSC_ERR_ARG_WRONG, "Non-order object");

    PetscHeaderCreate(o, _VarOrdering, int, VAR_ORDER_COOKIE, 0, "VarOrdering", map->comm, VarOrderingDestroy, 0);
    PetscLogObjectCreate(o);
    ierr = PetscObjectCompose((PetscObject) o, "ClassMap", (PetscObject) map);                            CHKERRQ(ierr);

    ierr = PetscBinaryRead(fd, &o->numVars,            1,               PETSC_INT);                       CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &o->numLocVars,         1,               PETSC_INT);                       CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &o->numOverlapVars,     1,               PETSC_INT);                       CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &o->numNewVars,         1,               PETSC_INT);                       CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &o->numLocNewVars,      1,               PETSC_INT);                       CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &o->numTotalFields,     1,               PETSC_INT);                       CHKERRQ(ierr);
    ierr = PetscMalloc((numProcs+1)      * sizeof(int),   &o->firstVar);                                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd,  o->firstVar,           numProcs+1,      PETSC_INT);                       CHKERRQ(ierr);
    ierr = PetscMalloc(numOverlapNodes   * sizeof(int),   &o->offsets);                                   CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd,  o->offsets,            numOverlapNodes, PETSC_INT);                       CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &hasLocOffsets,         1,               PETSC_INT);                       CHKERRQ(ierr);
    if (hasLocOffsets) {
      ierr = PetscMalloc(numGhostNodes   * sizeof(int),   &o->localOffsets);                              CHKERRQ(ierr);
      ierr = PetscBinaryRead(fd,  o->localOffsets,     numGhostNodes,   PETSC_INT);                       CHKERRQ(ierr);
    }
    ierr = PetscMalloc(o->numTotalFields * sizeof(int *), &o->localStart);                                CHKERRQ(ierr);
    for(f = 0; f < map->numFields; f++) {
      field = map->fields[f];
      ierr = PetscMalloc(numClasses      * sizeof(int),   &o->localStart[field]);                         CHKERRQ(ierr);
      ierr = PetscBinaryRead(fd, o->localStart[field], numClasses,      PETSC_INT);                       CHKERRQ(ierr);
    }
    PetscLogObjectMemory(o, (numProcs+1 + numOverlapNodes + o->numTotalFields*numClasses) * sizeof(int)
                     + o->numTotalFields * sizeof(int *));
    *order = o;
  }

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingDuplicate"
/*@C VarOrderingDuplicate
  This function duplicates a global variable ordering.

  Collective on VarOrdering

  Input Parameter:
. order    - The ordering

  Output Parameter:
. neworder - The new ordering

  Level: beginner

.keywords variable ordering, finite element
.seealso VarOrderingDuplicateIndices()
@*/
int VarOrderingDuplicate(VarOrdering order, VarOrdering *neworder)
{
  VarOrdering   o;
  FieldClassMap map;
  int           numFields;
  int          *fields;
  int           numOverlapNodes;
  int           numGhostNodes;
  int           numClasses;
  int           numProcs;
  int           field, newField;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  PetscValidPointer(neworder);

  ierr = MPI_Comm_size(order->comm, &numProcs);                                                           CHKERRQ(ierr);

  /* Create the ordering */
  PetscHeaderCreate(o, _VarOrdering, int, VAR_ORDER_COOKIE, 0, "VarOrdering", order->comm, VarOrderingDestroy, 0);
  PetscLogObjectCreate(o);
  ierr = PetscObjectQuery((PetscObject) order, "ClassMap", (PetscObject *) &map);                         CHKERRQ(ierr);
  ierr = PetscObjectCompose((PetscObject) o, "ClassMap", (PetscObject) map);                              CHKERRQ(ierr);
  numFields       = map->numFields;
  fields          = map->fields;
  numOverlapNodes = map->numOverlapNodes;
  numGhostNodes   = map->numGhostNodes;
  numClasses      = map->numClasses;

  /* Set sizes */
  o->numVars           = order->numVars;
  o->numLocVars        = order->numLocVars;
  o->numOverlapVars    = order->numOverlapVars;
  o->numNewVars        = order->numNewVars;
  o->numLocNewVars     = order->numLocNewVars;
  o->numOverlapNewVars = order->numOverlapNewVars;
  o->numTotalFields    = order->numTotalFields;
  /* Allocate memory */
  ierr = PetscMalloc((numProcs+1)    * sizeof(int), &o->firstVar);                                        CHKERRQ(ierr);
  ierr = PetscMalloc(numOverlapNodes * sizeof(int), &o->offsets);                                         CHKERRQ(ierr);
  /* Copy data */
  ierr = PetscMemcpy(o->firstVar, order->firstVar, (numProcs+1)    * sizeof(int));                        CHKERRQ(ierr);
  ierr = PetscMemcpy(o->offsets,  order->offsets,  numOverlapNodes * sizeof(int));                        CHKERRQ(ierr);
  PetscLogObjectMemory(o, (numProcs+1 + numOverlapNodes) * sizeof(int));

  if (order->localOffsets != PETSC_NULL) {
    ierr = PetscMalloc(numGhostNodes * sizeof(int), &o->localOffsets);                                    CHKERRQ(ierr);
    ierr = PetscMemcpy(o->localOffsets, order->localOffsets, numGhostNodes * sizeof(int));                CHKERRQ(ierr);
    PetscLogObjectMemory(o, numGhostNodes * sizeof(int));
  }

  ierr = PetscMalloc(o->numTotalFields * sizeof(int *), &o->localStart);                                  CHKERRQ(ierr);
  for(field = 0; field < numFields; field++) {
    newField = fields[field];
    ierr = PetscMalloc(numClasses * sizeof(int), &o->localStart[newField]);                               CHKERRQ(ierr);
    ierr = PetscMemcpy(o->localStart[newField], order->localStart[newField], numClasses * sizeof(int));   CHKERRQ(ierr);
  }
  PetscLogObjectMemory(o, o->numTotalFields * sizeof(int *) + numClasses*numFields * sizeof(int));

  *neworder = o;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingDuplicateIndices"
/*@C VarOrderingDuplicateIndices
  This function copies the global variable ordering indices to another ordering.

  Collective on VarOrdering

  Input Parameter:
. order    - The original ordering

  Output Parameter:
. neworder - The ordering with updated indices

  Level: intermediate

.keywords variable ordering, finite element
.seealso VarOrderingDuplicate()
@*/
int VarOrderingDuplicateIndices(VarOrdering order, VarOrdering neworder)
{
  PetscFunctionBegin;
  SETERRQ(PETSC_ERR_SUP, "Not yet implemented");
#if 0
  if ((mat->rowSize != newmat->rowSize) || (mat->colSize != newmat->colSize)) SETERRQ(1, "Incompatible matrix sizes");
  ierr = PetscMemcpy(newmat->rowIndices, mat->rowIndices, mat->rowSize * sizeof(int));                   CHKERRQ(ierr);
  ierr = PetscMemcpy(newmat->colIndices, mat->colIndices, mat->colSize * sizeof(int));                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
#endif
}

#undef  __FUNCT__
#define __FUNCT__ "VarOrderingCompatible"
/*@
  VarOrderingCompatible - This function determines whether the two orderings are shape compatible.

  Collective on VarOrdering

  Input Parameters:
+ order      - The variable ordering
- anOrder    - Another ordering

  Output Parameter:
. compatible - The compatibility flag

  Level: intermediate

.keywords: variable ordering, compatible
.seealso: VarOrderingGetLocalSize()
@*/
int VarOrderingCompatible(VarOrdering order, VarOrdering anOrder, PetscTruth *compatible)
{
  int locCompat = 1;
  int compat;
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(order,   VAR_ORDER_COOKIE);
  PetscValidHeaderSpecific(anOrder, VAR_ORDER_COOKIE);
  PetscValidPointer(compatible);
  if ((order->numVars        != anOrder->numVars)    ||
      (order->numLocVars     != anOrder->numLocVars) ||
      (order->numOverlapVars != anOrder->numOverlapVars)) {
    locCompat = 0;
  }
  ierr = MPI_Allreduce(&locCompat, &compat, 1, MPI_INT, MPI_LAND, order->comm);                           CHKERRQ(ierr);
  if (compat) {
    *compatible = PETSC_TRUE;
  } else {
    *compatible = PETSC_FALSE;
  }
  PetscFunctionReturn(0);
}

/*---------------------------------------------- Local Variable Ordering ---------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingCreate"
/*@C LocalVarOrderingCreate
  This function creates a local variable ordering for a finite element calculation.

  Collective on Grid

  Input Parameter:
+ grid      - The underlying grid for the ordering
. numFields - The number of fields in the ordering
- fields    - The fields in the ordering

  Output Parameter :
. order     - The ordering

  Level: beginner

.keywords variable ordering, finite element
.seealso LocalVarOrderingDestroy
@*/
int LocalVarOrderingCreate(Grid grid, int numFields, int *fields, LocalVarOrdering *order)
{ 
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidIntPointer(fields);
  PetscValidPointer(order);
  ierr = GridSetUp(grid);                                                                                CHKERRQ(ierr);
  ierr = (*grid->ops->gridcreatelocalvarordering)(grid, numFields, fields, order);                       CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingDestroy"
/*@C LocalVarOrderingDestroy
  This function destroys a local variable ordering.

  Not collective

  Input Parameter:
. order - The ordering

  Level: beginner

.keywords variable ordering, finite element
.seealso LocalVarOrderingCreate
@*/
int LocalVarOrderingDestroy(LocalVarOrdering order)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(order, VAR_ORDER_COOKIE);
  ierr = PetscFree(order->elemStart);                                                                    CHKERRQ(ierr);
  ierr = PetscFree(order->fields);                                                                       CHKERRQ(ierr);
  PetscLogObjectDestroy(order);
  PetscHeaderDestroy(order);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingSerialize"
/*@ 
  LocalVarOrderingSerialize - This function stores or recreates a local variable ordering using a viewer for a binary file.

  Collective on Grid

  Input Parameters:
+ grid   - The grid for the ordering
. viewer - The viewer context
- store  - This flag is PETSC_TRUE is data is being written, otherwise it will be read

  Output Parameter:
. order  - The ordering

  Level: beginner

.keywords: grid vector, serialize
.seealso: GridSerialize()
@*/
int LocalVarOrderingSerialize(Grid grid, LocalVarOrdering *order, PetscViewer viewer, PetscTruth store)
{
  LocalVarOrdering o;
  MPI_Comm         comm;
  int              fd;
  int              cookie, numFields;
  PetscTruth       match;
  int              ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE);
  PetscValidPointer(order);

  ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &match);                             CHKERRQ(ierr);
  if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONG, "Must be binary viewer");
  ierr = PetscViewerBinaryGetDescriptor(viewer, &fd);                                                     CHKERRQ(ierr);
  ierr = GridGetNumFields(grid, &numFields);                                                              CHKERRQ(ierr);

  if (store) {
    o = *order;
    PetscValidHeaderSpecific(o, VAR_ORDER_COOKIE);
    ierr = PetscBinaryWrite(fd, &o->cookie,    1,             PETSC_INT,  0);                             CHKERRQ(ierr);

    ierr = PetscBinaryWrite(fd, &o->numFields, 1,             PETSC_INT,  0);                             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  o->fields,    o->numFields,  PETSC_INT,  0);                             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &o->elemSize,  1,             PETSC_INT,  0);                             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  o->elemStart, numFields,     PETSC_INT,  0);                             CHKERRQ(ierr);
  } else {
    ierr = PetscBinaryRead(fd, &cookie,       1,             PETSC_INT);                                  CHKERRQ(ierr);
    if (cookie != VAR_ORDER_COOKIE) SETERRQ(PETSC_ERR_ARG_WRONG, "Non-order object");

    ierr = PetscObjectGetComm((PetscObject) grid, &comm);                                                 CHKERRQ(ierr);
    PetscHeaderCreate(o, _LocalVarOrdering, int, VAR_ORDER_COOKIE, 0, "VarOrdering", comm, LocalVarOrderingDestroy, PETSC_NULL);
    PetscLogObjectCreate(o);
    ierr = PetscBinaryRead(fd, &o->numFields, 1,             PETSC_INT);                                  CHKERRQ(ierr);

    ierr = PetscMalloc(o->numFields * sizeof(int), &o->fields);                                           CHKERRQ(ierr);
    ierr = PetscMalloc(numFields    * sizeof(int), &o->elemStart);                                        CHKERRQ(ierr);
    PetscLogObjectMemory(o, (o->numFields + numFields) * sizeof(int));
    ierr = PetscBinaryRead(fd,  o->fields,    o->numFields,  PETSC_INT);                                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &o->elemSize,  1,             PETSC_INT);                                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd,  o->elemStart, numFields,     PETSC_INT);                                  CHKERRQ(ierr);
    *order = o;
  }

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingDuplicate"
/*@C LocalVarOrderingDuplicate
  This function duplicates a local variable ordering.

  Collective on LocalVarOrdering

  Input Parameter:
. order    - The ordering

  Output Parameter:
. neworder - The new ordering

  Level: beginner

.keywords variable ordering, finite element
.seealso LocalVarOrderingDuplicateIndices()
@*/
int LocalVarOrderingDuplicate(LocalVarOrdering order, LocalVarOrdering *neworder)
{
  SETERRQ(PETSC_ERR_SUP, "Not yet implemented");
#if 0
  PetscFunctionBegin;
  PetscFunctionReturn(0);
#endif
}

#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingDuplicateIndices"
/*@C LocalVarOrderingDuplicateIndices
  This function copies the local variable ordering indices to another ordering.

  Collective on LocalVarOrdering

  Input Parameter:
. order    - The original ordering

  Output Parameter:
. neworder - The ordering with updated indices

  Level: intermediate

.keywords variable ordering, finite element
.seealso LocalVarOrderingDuplicate()
@*/
int LocalVarOrderingDuplicateIndices(LocalVarOrdering order, LocalVarOrdering neworder)
{
  SETERRQ(PETSC_ERR_SUP, "Not yet implemented");
#if 0
  PetscFunctionBegin;
  if ((mat->rowSize != newmat->rowSize) || (mat->colSize != newmat->colSize)) SETERRQ(1, "Incompatible matrix sizes");
  ierr = PetscMemcpy(newmat->rowIndices, mat->rowIndices, mat->rowSize * sizeof(int));                   CHKERRQ(ierr);
  ierr = PetscMemcpy(newmat->colIndices, mat->colIndices, mat->colSize * sizeof(int));                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
#endif
}

#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingGetSize"
/*@C
  LocalVarOrderingGetSize - This function returns the greatest number of variables on any node,
  or equivalently the number of shape functions in the element matrix.

  Not collective

  Input Parameter:
. order - The ordering

  Output Parameter:
. size  - The size

  Level: intermediate

.keywords local, order, size
.seealso LocalVarOrderingGetFieldStart()
@*/
int LocalVarOrderingGetSize(LocalVarOrdering order, int *size) {
  PetscFunctionBegin;
  PetscValidIntPointer(size);
  *size = order->elemSize;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingGetField"
/*@C
  LocalVarOrderingGetField - This function returns the field at a given index in this ordering.

  Not collective

  Input Parameters:
+ order - The ordering
- f     - The field index

  Output Parameter:
. field - The field

  Level: intermediate

.keywords local, order, field
.seealso LocalVarOrderingGetFieldStart()
@*/
int LocalVarOrderingGetField(LocalVarOrdering order, int f, int *field) {
  PetscFunctionBegin;
  PetscValidIntPointer(field);
  if ((f < 0) || (f >= order->numFields)) {
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid field index %d must be in [0,%d)", f, order->numFields);
  }
  *field = order->fields[f];
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingGetFieldIndex"
/*@C
  LocalVarOrderingGetFieldIndex - This function returns the field index in this ordering of the given field.

  Not collective

  Input Parameters:
+ order - The ordering
- field - The field

  Output Parameter:
. f     - The field index

  Level: intermediate

.keywords local, order, field, index
.seealso LocalVarOrderingGetFieldStart()
@*/
int LocalVarOrderingGetFieldIndex(LocalVarOrdering order, int field, int *f) {
  int index;

  PetscFunctionBegin;
  PetscValidIntPointer(f);
  for(index = 0; index < order->numFields; index++) {
    if (order->fields[index] == field) break;
  }
  if (index >= order->numFields) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field %d not present in order", field);
  *f = index;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "LocalVarOrderingGetFieldStart"
/*@C
  LocalVarOrderingGetFieldStart - This function returns the first index of the given field in the element vector or matrix.

  Not collective

  Input Parameters:
+ order - The ordering
- field - The field

  Output Parameter:
. start - The element vector or matrix index, or -1 if the field is not active

  Level: intermediate

.keywords local, order, start, field
.seealso LocalVarGetSize()
@*/
int LocalVarOrderingGetFieldStart(LocalVarOrdering order, int field, int *start) {
  PetscFunctionBegin;
  PetscValidIntPointer(start);
  *start = order->elemStart[field];
  PetscFunctionReturn(0);
}