#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: gridDB.c,v 1.4 2000/07/16 05:40:27 knepley Exp $";
#endif

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

/*-------------------------------------------------- Field Functions ------------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "GridAddField"
/*@C
  GridAddField - This function defines a field on the grid.

  Not collective

  Input Parameters:
+ grid     - The grid
. name     - The field name
. discType - The DiscretizationType for the field
- numComp  - The nubmer of components in the field

  Output Parameter:
. field    - The canonical field number

  Level: intermediate

.keywords: grid, field,
.seealso: GridSetFieldName()
@*/
int GridAddField(Grid grid, const char name[], DiscretizationType discType, int numComp, int *field)
{
  Field *tempFields;
  int    ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  if (numComp <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of field components %d", numComp);
  while (grid->numFields >= grid->maxFields) {
    ierr = PetscMalloc(grid->maxFields*2 * sizeof(Field), &tempFields);                                   CHKERRQ(ierr);
    PetscLogObjectMemory(grid, grid->maxFields * sizeof(Field));
    ierr = PetscMemcpy(tempFields, grid->fields, grid->maxFields * sizeof(Field));                        CHKERRQ(ierr);
    ierr = PetscFree(grid->fields);                                                                       CHKERRQ(ierr);
    grid->fields     = tempFields;
    grid->maxFields *= 2;
  }
  if (field != PETSC_NULL) {
    PetscValidIntPointer(field);
    *field = grid->numFields;
  }

  ierr = PetscStrallocpy(name, &grid->fields[grid->numFields].name);                                      CHKERRQ(ierr);
  ierr = PetscStrallocpy(discType, &grid->fields[grid->numFields].discType);                              CHKERRQ(ierr);
  grid->fields[grid->numFields].numComp            = numComp;
  grid->fields[grid->numFields].isConstrained      = PETSC_FALSE;
  grid->fields[grid->numFields].constraintCompDiff = 0;
  grid->fields[grid->numFields].isActive           = PETSC_FALSE;
  ierr = DiscretizationCreate(grid->comm, &grid->fields[grid->numFields].disc);                           CHKERRQ(ierr);
  ierr = DiscretizationSetNumComponents(grid->fields[grid->numFields].disc, numComp);                     CHKERRQ(ierr);
  ierr = DiscretizationSetField(grid->fields[grid->numFields].disc, grid->numFields);                     CHKERRQ(ierr);
  ierr = DiscretizationSetType(grid->fields[grid->numFields].disc, discType);                             CHKERRQ(ierr);
  ierr = DiscretizationSetFromOptions(grid->fields[grid->numFields].disc);                                CHKERRQ(ierr);
  grid->numFields++;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetNumFields"
/*@C
  GridGetNumFields - This function gets the number of fields in the grid.

  Not collective

  Input Parameter:
. grid - The grid

  Output Parameter:
. num  - The number of fields

  Level: intermediate

.keywords: grid, field,
.seealso: GridSetFieldName()
@*/
int GridGetNumFields(Grid grid, int *num)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidIntPointer(num);
  *num = grid->numFields;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetFieldName"
/*@C
  GridGetFieldName - This function gets the name of a particular field.

  Not collective

  Input Parameters:
+ grid  - The grid
- field - The field to name

  Output Parameter:
. name  - The name string

  Level: intermediate

.keywords: grid, field, name
.seealso: GridSetFieldName()
@*/
int GridGetFieldName(Grid grid, int field, char **name)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidPointer(name);
  GridValidField(grid, field);
  if (grid->fields[field].name == PETSC_NULL) {
    *name = PETSC_NULL;
  } else {
    ierr = PetscStrallocpy(grid->fields[field].name, name);                                               CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridSetFieldName"
/*@C
  GridSetFieldName - This function sets the name of a particular field.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
. field - The field to name
- name  - The name string

  Level: intermediate

.keywords: grid, field, name
.seealso: GridGetFieldName()
@*/
int GridSetFieldName(Grid grid, int field, char *name)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, field);
  PetscValidCharPointer(name);
  ierr = PetscStrfree(grid->fields[field].name);                                                          CHKERRQ(ierr);
  ierr = PetscStrallocpy(name, &grid->fields[field].name);                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetFieldComponents"
/*@C
  GridGetFieldComponents - This function gets the number of components in a field.

  Not collective

  Input Parameters
+ grid  - The grid
- field - The field

  Output Parameter:
. comp - The number of components in the field

  Level: intermediate

.keywords: grid, field, component
.seealso: GridGetFieldDiscType()
@*/
int GridGetFieldComponents(Grid grid, int field, int *comp)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, field);
  PetscValidIntPointer(comp);
  *comp = grid->fields[field].numComp;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetFieldDisc"
/*@C
  GridGetFieldDisc - This function gets the discretization for a field.

  Not collective

  Input Parameters:
+ grid  - The grid
- field - The field

  Output Parameter:
. disc  - The discretization

  Level: intermediate

.keywords: grid, field, discretization
.seealso: GridGetFieldComponents()
@*/
int GridGetFieldDisc(Grid grid, int field, Discretization *disc)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidIntPointer(disc);
  GridValidField(grid, field);
  *disc = grid->fields[field].disc;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetNumActiveFields"
/*@C
  GridGetNumActiveFields - This function gets the number of fields participating in the calculation.

  Not collective

  Input Parameter:
. grid - The grid

  Output Parameter:
. num  - The number of active fields

  Level: intermediate

.keywords: grid, field, active
.seealso: GridSetFieldName()
@*/
int GridGetNumActiveFields(Grid grid, int *num)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidIntPointer(num);
  *num = grid->numActiveFields;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetActiveField"
/*@C
  GridGetActiveField - This function gets the canonical field number for a certain active field.

  Not collective

  Input Parameters:
+ grid  - The grid
- n     - The nth active field

  Output Parameter:
. field - The field

  Level: intermediate

.keywords: grid, field, active
.seealso: GridSetFieldName()
@*/
int GridGetActiveField(Grid grid, int n, int *field)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidIntPointer(field);
  if ((n < 0) || (n >= grid->numActiveFields)) {
    SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid active field number %d should be in [0,%d)", n, grid->numActiveFields);
  }
  *field = grid->defaultFields[n];
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridSetActiveField"
/*@C
  GridSetActiveField - This function makes one field active for a particular grid.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
- field - The field to make active

  Level: intermediate

.keywords: grid, active, field
.seealso: GridAddActiveField()
@*/
int GridSetActiveField(Grid grid, int field)
{
  int fieldIndex;
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, field);
  grid->numActiveFields = 0;
  for(fieldIndex = 0; fieldIndex < grid->numFields; fieldIndex++) grid->fields[fieldIndex].isActive = PETSC_FALSE;
  ierr = GridAddActiveField(grid, field);                                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridAddActiveField"
/*@C
  GridAddActiveField - This function makes another field active for a particular grid.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
- field - The field to make active

  Level: intermediate

.keywords: grid, active, field
.seealso: GridSetActiveField()
@*/
int GridAddActiveField(Grid grid, int field)
{
  int *tempFields;
  int  ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE); 
  GridValidField(grid, field);
  /* Setup assembly variables */
  while (grid->numActiveFields >= grid->maxActiveFields) {
    ierr = PetscMalloc(grid->maxActiveFields*2 * sizeof(int), &tempFields);                               CHKERRQ(ierr);
    PetscLogObjectMemory(grid, grid->maxActiveFields * sizeof(int));
    ierr = PetscMemcpy(tempFields, grid->defaultFields, grid->maxActiveFields * sizeof(int));             CHKERRQ(ierr);
    ierr = PetscFree(grid->defaultFields);                                                                CHKERRQ(ierr);
    grid->defaultFields = tempFields;
    grid->maxActiveFields *= 2;
  }
  if (grid->fields[field].isActive == PETSC_FALSE) grid->defaultFields[grid->numActiveFields++] = field;
  grid->fields[field].isActive = PETSC_TRUE;
  grid->setupcalled            = PETSC_FALSE;
  grid->bdSetupCalled          = PETSC_FALSE;
  PetscFunctionReturn(0);
}

/*------------------------------------------------ Variable Functions -----------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "GridGetNumVars"
/*@C GridGetNumVars
  This function gets the number of variables in the grid.

  Not collective

  Input Parameter:
. grid - The grid

  Output Parameters:
+ locNum      - The number of local variables
. num         - The number of variables
. locConstNum - The number of local constrained variables
- constNum    - The number of constrained variables

  Level: intermediate

.keywords grid, field,
.seealso GridSetFieldName
@*/
int GridGetNumVars(Grid grid, int *locNum, int *num, int *locConstNum, int *constNum)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  if (grid->order == PETSC_NULL) {
    ierr = GridSetUp(grid);                                                                              CHKERRQ(ierr);
  }
  if (locNum != PETSC_NULL) {
    PetscValidIntPointer(locNum);
    *locNum = grid->order->numLocVars;
  }
  if (num    != PETSC_NULL) {
    PetscValidIntPointer(num);
    *num    = grid->order->numVars;
  }
  if (grid->isConstrained == PETSC_TRUE) {
    if ((grid->constraintOrder == PETSC_NULL) && ((locConstNum != PETSC_NULL) || (constNum != PETSC_NULL))) {
      SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Constraints are not yet setup");
    }
    if (locConstNum != PETSC_NULL) {
      PetscValidIntPointer(locConstNum);
      *locConstNum = grid->constraintOrder->numLocVars;
    }
    if (constNum    != PETSC_NULL) {
      PetscValidIntPointer(constNum);
      *constNum    = grid->constraintOrder->numVars;
    }
  }
  PetscFunctionReturn(0);
}

/*-------------------------------------------------- Rhs Functions --------------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "GridGetNumNonlinearOperators"
/*@
  GridGetNumNonlinearOperators - Returns the number of nonlinear operators for the problem defined on the grid.

  Not collective

  Input Parameter:
. grid   - The grid

  Output Parameter:
. numOps - The number of operators

  Level: intermediate

.keywords: grid, operator, nonlinear
.seealso: GridGetNonlinearOperatorStart()
@*/
int GridGetNumNonlinearOperators(Grid grid, int *numOps)
{
  int op;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidIntPointer(numOps);
  *numOps = 0;
  for(op = 0; op < grid->numRhsOps; op++) {
    if (grid->rhsOps[op].nonlinearOp != PETSC_NULL) *numOps += 1;
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetNonlinearOperatorStart"
/*@
  GridGetNonlinearOperatorStart - Retrieves information about the first nonlinear
  operator for the problem defined on the grid.

  Not collective

  Input Parameter:
. grid   - The grid

  Output Parameters:
+ op     - The operator
. field  - The shape function field
. alpha  - The multiplier of this operator
- isALE  - The flag for ALE operators

  Level: intermediate

.keywords: grid, operator, nonlinear
.seealso: GridGetNonlinearOperatorNext(), GridGetNumNonlinearOperators()
@*/
int GridGetNonlinearOperatorStart(Grid grid, NonlinearOperator *op, int *field, PetscScalar *alpha, PetscTruth *isALE)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  grid->activeNonlinearOp = -1;
  ierr = GridGetNonlinearOperatorNext(grid, op, field, alpha, isALE);                                    CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetNonlinearOperatorNext"
/*@
  GridGetNonlinearOperatorNext - Retrieves information about the next nonlinear
  operator for the problem defined on the grid.

  Not collective

  Input Parameter:
. grid   - The grid

  Output Parameters:
+ op     - The operator
. field  - The shape function field
. alpha  - The multiplier of this operator
- isALE  - The flag for ALE operators

  Level: intermediate

.keywords: grid, operator, nonlinear
.seealso: GridNonlinearGetOperatorStart(), GridGetNumNonlinearOperators()
@*/
int GridGetNonlinearOperatorNext(Grid grid, NonlinearOperator *op, int *field, PetscScalar *alpha, PetscTruth *isALE)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  grid->activeNonlinearOp++;
  while((grid->rhsOps[grid->activeNonlinearOp].nonlinearOp == PETSC_NULL) && (grid->activeNonlinearOp < grid->numRhsOps)) {
    grid->activeNonlinearOp++;
  }
  if (grid->activeNonlinearOp >= grid->numRhsOps) {
    if (op    != PETSC_NULL) *op    = PETSC_NULL;
    if (field != PETSC_NULL) *field = -1;
    if (alpha != PETSC_NULL) *alpha = 0.0;
    if (isALE != PETSC_NULL) *isALE = PETSC_FALSE;
    grid->activeNonlinearOp = -1;
    PetscFunctionReturn(0);
  }
  if (op    != PETSC_NULL) {
    PetscValidIntPointer(op);
    *op    = grid->rhsOps[grid->activeNonlinearOp].nonlinearOp;
  }
  if (field != PETSC_NULL) {
    PetscValidIntPointer(field);
    *field = grid->rhsOps[grid->activeNonlinearOp].field;
  }
  if (alpha != PETSC_NULL) {
    PetscValidScalarPointer(alpha);
    *alpha = grid->rhsOps[grid->activeNonlinearOp].alpha;
  }
  if (isALE != PETSC_NULL) {
    PetscValidIntPointer(isALE);
    *isALE = grid->rhsOps[grid->activeNonlinearOp].isALE;
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridAddRhsFunction"
/*@C GridAddRhsFunction
  This function adds the point function to use on the rhs for the given field.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
. field - The field
. func  - The point function
- alpha - A scalar multiplier

  Level: beginner

.keywords active field
.seealso GridSetMatOperator, GridSetRhsOperator
@*/
int GridAddRhsFunction(Grid grid, int field, PointFunction f, PetscScalar alpha)
{
  GridFunc *tempFuncs;
  int       ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, field);
  while (grid->numRhsFuncs >= grid->maxRhsFuncs) {
    ierr = PetscMalloc(grid->maxRhsFuncs*2 * sizeof(GridOp), &tempFuncs);                                 CHKERRQ(ierr);
    PetscLogObjectMemory(grid, grid->maxRhsFuncs * sizeof(GridOp));
    ierr = PetscMemcpy(tempFuncs, grid->rhsFuncs, grid->maxRhsFuncs * sizeof(GridOp));                    CHKERRQ(ierr);
    ierr = PetscFree(grid->rhsFuncs);                                                                     CHKERRQ(ierr);
    grid->rhsFuncs     = tempFuncs;
    grid->maxRhsFuncs *= 2;
  }
  grid->rhsFuncs[grid->numRhsFuncs].func  = f;
  grid->rhsFuncs[grid->numRhsFuncs].field = field;
  grid->rhsFuncs[grid->numRhsFuncs].alpha = alpha;
  grid->numRhsFuncs++;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridAddRhsNonlinearOperator"
/*@C GridAddRhsNonlinearOperator
  This function add the nonlinear operator to use on the Rhs for the given field.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
. field - The field
. f     - The function defining the nonlinear operator
. alpha - A scalar multiplier
- isALE - A flag for ALE operators

  Level: beginner

.keywords active field
.seealso GridSetRhsOperator, GridAddRhsOperator
@*/
int GridAddRhsNonlinearOperator(Grid grid, int field, NonlinearOperator f, PetscScalar alpha, PetscTruth isALE)
{
  GridOp *tempOps;
  int     ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, field);
  while (grid->numRhsOps >= grid->maxRhsOps) {
    ierr = PetscMalloc(grid->maxRhsOps*2 * sizeof(GridOp), &tempOps);                                     CHKERRQ(ierr);
    PetscLogObjectMemory(grid, grid->maxRhsOps * sizeof(GridOp));
    ierr = PetscMemcpy(tempOps, grid->rhsOps, grid->maxRhsOps * sizeof(GridOp));                          CHKERRQ(ierr);
    ierr = PetscFree(grid->rhsOps);                                                                       CHKERRQ(ierr);
    grid->rhsOps     = tempOps;
    grid->maxRhsOps *= 2;
  }
  grid->rhsOps[grid->numRhsOps].nonlinearOp = f;
  grid->rhsOps[grid->numRhsOps].field       = field;
  grid->rhsOps[grid->numRhsOps].alpha       = alpha;
  grid->rhsOps[grid->numRhsOps].isALE       = isALE;
  grid->rhsOps[grid->numRhsOps].op          = -1;
  grid->numRhsOps++;
  if (isALE == PETSC_TRUE) grid->ALEActive = PETSC_TRUE;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridSetRhsOperator"
/*@C GridSetRhsOperator
  This function sets the operator to use on the Rhs for the given field.

  Collective on Grid

  Input Parameters:
+ grid   - The grid
. sfield - The field for shape functions
. tfield - The field for test functions
. op     - The operator
. alpha  - A scalar multiple for the operator
- isALE  - A flag for ALE operators

  Output Parameter:
. index  - [Optional] A unique indentifier for the operator

  Level: beginner

.keywords active field
.seealso GridAddRhsOperator, GridAddMatOperator, GridSetRhsFunction
@*/
int GridSetRhsOperator(Grid grid, int sField, int tField, int op, PetscScalar alpha, PetscTruth isALE, int *index)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  grid->numRhsOps = 0;
  ierr = GridAddRhsOperator(grid, sField, tField, op, alpha, isALE, index);                              CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridAddRhsOperator"
/*@C GridAddRhsOperator
  This function adds the operator to the Rhs for the given field.

  Collective on Grid

  Input Parameters:
+ grid   - The grid
. sfield - The field for shape functions
. tfield - The field for test functions
. op     - The operator
. alpha  - A scalar multiple for the operator
- isALE  - A flag for ALE operators

  Output Parameter:
. index  - [Optional] A unique indentifier for the operator

  Level: beginner

.keywords active field
.seealso GridSetRhsOperator, GridSetMatOperator, GridSetRhsFunction
@*/
int GridAddRhsOperator(Grid grid, int sField, int tField, int op, PetscScalar alpha, PetscTruth isALE, int *index)
{
  Discretization sDisc = grid->fields[sField].disc;
  Discretization tDisc = grid->fields[tField].disc;
  int            dim   = grid->dim;
  GridOp        *tempOps;
  int            ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, sField);
  GridValidField(grid, tField);
  if ((op < 0) || (op >= sDisc->numOps) || (!sDisc->operators[op])) {
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d", op);
  }
  if ((sDisc->operators[op]->precompInt == PETSC_NULL) &&
      (sDisc->operators[op]->opFunc     == PETSC_NULL) &&
      (sDisc->operators[op]->ALEOpFunc  == PETSC_NULL)) {
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d has no associated function", op);
  }
  if (sDisc->numQuadPoints != tDisc->numQuadPoints) {
    SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets");
  }
  while (grid->numRhsOps >= grid->maxRhsOps) {
    ierr = PetscMalloc(grid->maxRhsOps*2 * sizeof(GridOp), &tempOps);                                     CHKERRQ(ierr);
    PetscLogObjectMemory(grid, grid->maxRhsOps * sizeof(GridOp));
    ierr = PetscMemcpy(tempOps, grid->rhsOps, grid->maxRhsOps * sizeof(GridOp));                          CHKERRQ(ierr);
    ierr = PetscFree(grid->rhsOps);                                                                       CHKERRQ(ierr);
    grid->rhsOps     = tempOps;
    grid->maxRhsOps *= 2;
  }
  if (index != PETSC_NULL) {
    PetscValidIntPointer(index);
    *index = grid->numRhsOps;
  }
  grid->rhsOps[grid->numRhsOps].op          = op;
  grid->rhsOps[grid->numRhsOps].field       = sField;
  grid->rhsOps[grid->numRhsOps].alpha       = alpha;
  grid->rhsOps[grid->numRhsOps].isALE       = isALE;
  grid->rhsOps[grid->numRhsOps].nonlinearOp = PETSC_NULL;
  grid->numRhsOps++;
  if (isALE == PETSC_TRUE) grid->ALEActive = PETSC_TRUE;
  /* Assign test function field to operator -- Need more sophisticated checking */
  switch (op) {
  case IDENTITY:
  case LAPLACIAN:
    if (sDisc->size != tDisc->size) {
      SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
    }
    break;
  case GRADIENT:
    if (sDisc->comp*dim != tDisc->comp) {
      SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
    }
    break;
  case DIVERGENCE:
    if (tDisc->comp*dim != sDisc->comp) {
      SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
    }
    break;
  }
  sDisc->operators[op]->test = tDisc;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridScaleRhs"
/*@C GridScaleRhs
  This function multiplies all the current components of the Rhs by the same scalar.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
- alpha - The scalar

  Level: intermediate

.keywords grid, rhs, scale
.seealso GridScaleSystemMatrix()
@*/
int GridScaleRhs(Grid grid, PetscScalar alpha)
{
  int func, op;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  /* Scale functions */
  for(func = 0; func < grid->numRhsFuncs; func++) grid->rhsFuncs[func].alpha *= alpha;
  /* Scale linear operators */
  for(op = 0; op < grid->numRhsOps; op++) grid->rhsOps[op].alpha *= alpha;
  PetscLogFlops(grid->numRhsFuncs + grid->numRhsOps);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridScaleRhsOperator"
/*@C GridScaleRhsOperator
  This function multiplies an operator of the Rhs by the a scalar.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
. alpha - The scalar
- index - The unique operator index returned by GridAddRhsOperator()

  Level: intermediate

.keywords grid, rhs, scale, operator
.seealso GridScaleSystemMatrix(), GridAddRhsOperator()
@*/
int GridScaleRhsOperator(Grid grid, PetscScalar alpha, int index)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  grid->rhsOps[index].alpha *= alpha;
  PetscLogFlops(1);
  PetscFunctionReturn(0);
}

/*--------------------------------------------- System Matrix Functions ---------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "GridIsMatrixFree"
/*@
  GridIsMatrixFree - Signals whether operators are applied matrix-free.

  Not collective

  Input Parameter:
. grid         - The grid

  Output Parameter:
. isMatrixFree - The flag for matrix-free operators

  Level: intermediate

.keywords: grid, operator, matrix-free
.seealso: GridSetMatrixFree()
@*/
int GridIsMatrixFree(Grid grid, PetscTruth *isMatrixFree)
{
  PetscFunctionBegin;
  *isMatrixFree = grid->isMatrixFree;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridSetMatrixFree"
/*@
  GridSetMatrixFree - Determines whether operators are applied matrix-free.

  Not collective

  Input Parameters:
. grid         - The grid
. isMatrixFree - The flag for matrix-free operators

  Level: intermediate

.keywords: grid, operator, matrix-free
.seealso: GridIsMatrixFree()
@*/
int GridSetMatrixFree(Grid grid, PetscTruth isMatrixFree)
{
  PetscFunctionBegin;
  grid->isMatrixFree = isMatrixFree;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetNumRhsOperators"
/*@
  GridGetNumRhsOperators - Returns the number of rhs operators for the problem defined on the grid.

  Not collective

  Input Parameter:
. grid   - The grid

  Output Parameter:
. numOps - The number of operators

  Level: intermediate

.keywords: grid, operator
.seealso: GridGetRhsOperatorStart()
@*/
int GridGetNumRhsOperators(Grid grid, int *numOps)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidIntPointer(numOps);
  *numOps = grid->numRhsOps;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetNumMatOperators"
/*@
  GridGetNumMatOperators - Returns the number of matrix operators for the problem defined on the grid.

  Not collective

  Input Parameter:
. grid   - The grid

  Output Parameter:
. numOps - The number of operators

  Level: intermediate

.keywords: grid, operator
.seealso: GridGetMatOperatorStart()
@*/
int GridGetNumMatOperators(Grid grid, int *numOps)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidIntPointer(numOps);
  *numOps = grid->numMatOps;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetMatOperatorStart"
/*@
  GridGetMatOperatorStart - Retrieves information about the first matrix operator
  for the problem defined on the grid.

  Not collective

  Input Parameter:
. grid   - The grid

  Output Parameters:
+ op     - The operator
. sField - The shape function field
. tField - The test function field
. alpha  - The multiplier of this operator
- isALE  - The flag for ALE operators

  Level: intermediate

.keywords: grid, operator
.seealso: GridGetMatOperatorNext(), GridGetNumMatOperators()
@*/
int GridGetMatOperatorStart(Grid grid, int *op, int *sField, int *tField, PetscScalar *alpha, PetscTruth *isALE)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  grid->activeMatOp = -1;
  ierr = GridGetMatOperatorNext(grid, op, sField, tField, alpha, isALE);                                 CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridGetMatOperatorNext"
/*@
  GridGetMatOperatorNext - Retrieves information about the next matrix operator
  for the problem defined on the grid.

  Not collective

  Input Parameter:
. grid   - The grid

  Output Parameters:
+ op     - The operator
. sField - The shape function field
. tField - The test function field
. alpha  - The multiplier of this operator
- isALE  - The flag for ALE operators

  Level: intermediate

.keywords: grid, operator
.seealso: GridMatGetOperatorStart(), GridGetNumMatOperators()
@*/
int GridGetMatOperatorNext(Grid grid, int *op, int *sField, int *tField, PetscScalar *alpha, PetscTruth *isALE)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  grid->activeMatOp++;
  if (grid->activeMatOp >= grid->numMatOps) {
    if (op     != PETSC_NULL) *op     = -1;
    if (sField != PETSC_NULL) *sField = -1;
    if (tField != PETSC_NULL) *tField = -1;
    if (alpha  != PETSC_NULL) *alpha  = 0.0;
    if (isALE  != PETSC_NULL) *isALE  = PETSC_FALSE;
    grid->activeMatOp = -1;
    PetscFunctionReturn(0);
  }
  if (op     != PETSC_NULL) {
    PetscValidIntPointer(op);
    *op     = grid->matOps[grid->activeMatOp].op;
  }
  if (sField != PETSC_NULL) {
    PetscValidIntPointer(sField);
    *sField = grid->matOps[grid->activeMatOp].field;
  }
  if (tField != PETSC_NULL) {
    PetscValidIntPointer(tField);
    *tField = grid->fields[grid->matOps[grid->activeMatOp].field].disc->operators[grid->matOps[grid->activeMatOp].op]->test->field;
  }
  if (alpha  != PETSC_NULL) {
    PetscValidScalarPointer(alpha);
    *alpha  = grid->matOps[grid->activeMatOp].alpha;
  }
  if (isALE  != PETSC_NULL) {
    PetscValidIntPointer(isALE);
    *isALE  = grid->matOps[grid->activeMatOp].isALE;
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridSetMatOperator"
/*@C GridSetMatOperator
  This function sets the operator to use in the system matrix for the given field.

  Collective on Grid

  Input Parameters:
+ grid   - The grid
. sField - The field for shape functions
. tField - The field for test functions
. op     - The operator
. alpha  - A scalar multiple for the operator
- isALE  - A flag for ALE operators

  Output Parameter:
. index  - [Optional] A unique indentifier for the operator

  Level: beginner

.keywords active field
.seealso GridAddMatOperator(), GridRemoveMatOperator(), GridAddRhsOperator()
@*/
int GridSetMatOperator(Grid grid, int sField, int tField, int op, PetscScalar alpha, PetscTruth isALE, int *index)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  grid->numMatOps = 0;
  ierr = GridAddMatOperator(grid, sField, tField, op, alpha, isALE, index);                              CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridAddMatOperator"
/*@C GridAddMatOperator
  This function adds the operator to the system matrix for the given field.

  Collective on Grid

  Input Parameters:
+ grid   - The grid
. sField - The field for shape functions
. tField - The field for test functions
. op     - The operator
. alpha  - A scalar multiple for the operator
- isALE  - A flag for ALE operators

  Output Parameter:
. index  - [Optional] A unique indentifier for the operator

  Level: beginner

.keywords active field
.seealso GridRemoveMatOperator(), GridSetMatOperator(), GridSetRhsOperator()
@*/
int GridAddMatOperator(Grid grid, int sField, int tField, int op, PetscScalar alpha, PetscTruth isALE, int *index)
{
  Discretization sDisc = grid->fields[sField].disc;
  Discretization tDisc = grid->fields[tField].disc;
  int            dim   = grid->dim;
  GridOp        *tempOps;
  int            ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, sField);
  GridValidField(grid, tField);
  if ((op < 0) || (op >= sDisc->numOps) || (!sDisc->operators[op])) {
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d", op);
  }
  if ((sDisc->operators[op]->precompInt == PETSC_NULL) &&
      (sDisc->operators[op]->opFunc     == PETSC_NULL) &&
      (sDisc->operators[op]->ALEOpFunc  == PETSC_NULL)) {
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d has no associated function", op);
  }
  if (sDisc->numQuadPoints != tDisc->numQuadPoints) {
    SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets");
  }
  while (grid->numMatOps >= grid->maxMatOps) {
    ierr = PetscMalloc(grid->maxMatOps*2 * sizeof(GridOp), &tempOps);                                     CHKERRQ(ierr);
    PetscLogObjectMemory(grid, grid->maxMatOps * sizeof(GridOp));
    ierr = PetscMemcpy(tempOps, grid->matOps, grid->maxMatOps * sizeof(GridOp));                          CHKERRQ(ierr);
    ierr = PetscFree(grid->matOps);                                                                       CHKERRQ(ierr);
    grid->matOps     = tempOps;
    grid->maxMatOps *= 2;
  }
  if (index != PETSC_NULL) {
    PetscValidIntPointer(index);
    *index = grid->numMatOps;
  }
  grid->matOps[grid->numMatOps].op          = op;
  grid->matOps[grid->numMatOps].field       = sField;
  grid->matOps[grid->numMatOps].alpha       = alpha;
  grid->matOps[grid->numMatOps].isALE       = isALE;
  grid->matOps[grid->numMatOps].nonlinearOp = PETSC_NULL;
  grid->numMatOps++;
  if (isALE == PETSC_TRUE) grid->ALEActive = PETSC_TRUE;
  /* Assign test function field to operator -- Need more sophisticated checking */
  switch (op) {
  case IDENTITY:
  case LAPLACIAN:
    if (sDisc->size != tDisc->size) {
      SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
    }
    break;
  case GRADIENT:
    if (sDisc->comp*dim != tDisc->comp) {
      SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
    }
    break;
  case DIVERGENCE:
    if (tDisc->comp*dim != sDisc->comp) {
      SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function field");
    }
    break;
  }
  sDisc->operators[op]->test = tDisc;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridRemoveMatOperator"
/*@C GridRemoveMatOperator
  This function removes the operator to the system matrix for the given field.

  Collective on Grid

  Input Parameters:
+ grid   - The grid
- index  - The operator ID from GridAddMatOperator()

  Level: beginner

.keywords active field
.seealso GridAddMatOperator(), GridSetMatOperator(), GridAddRhsOperator()
@*/
int GridRemoveMatOperator(Grid grid, int index) {
  int o;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  if ((index < 0) || (index >= grid->numMatOps)) {
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid operator ID %d should be in [0, %d)", index, grid->numMatOps);
  }
  for(o = index; o < grid->numMatOps-1; o++) {
    grid->matOps[o] = grid->matOps[o+1];
  }
  grid->numMatOps--;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridRegisterOperator"
/*@C GridRegisterOperator
  This function defines a new operator for a certain field

  Collective on Grid

  Input Parameters:
+ grid  - The grid
. field - The field on which the operator acts
- func  - The function deinfing the operator

  Output Parameter:
. newOp - The index of the new operator

  Level: intermediate

.keywords grid, operator, register
.seealso GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
@*/
int GridRegisterOperator(Grid grid, int field, OperatorFunction func, int *newOp)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, field);
  ierr = DiscretizationRegisterOperator(grid->fields[field].disc, func, newOp);                           CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridRegisterALEOperator"
/*@C GridRegisterALEOperator
  This function defines a new operator for a certain field

  Collective on Grid

  Input Parameters:
+ grid  - The grid
. field - The field on which the operator acts
- func  - The function deinfing the operator

  Output Parameter:
. newOp - The index of the new operator

  Level: intermediate

.keywords grid, oeprator, register
.seealso GridAddMatOperator, GridAddRhsOperator, GridSetRhsFunction
@*/
int GridRegisterALEOperator(Grid grid, int field, ALEOperatorFunction func, int *newOp)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  GridValidField(grid, field);
  ierr = DiscretizationRegisterALEOperator(grid->fields[field].disc, func, newOp);                        CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridScaleSystemMatrix"
/*@C GridScaleSystemMatrix
  This function multiplies all the current components of the system matrix by the same scalar.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
- alpha - The scalar

  Level: intermediate

.keywords grid, matrix, scale
.seealso GridScaleRhs()
@*/
int GridScaleSystemMatrix(Grid grid, PetscScalar alpha)
{
  int op;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  /* Scale linear operators */
  for(op = 0; op < grid->numMatOps; op++) grid->matOps[op].alpha *= alpha;
  PetscLogFlops(grid->numMatOps);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridScaleMatOperator"
/*@C GridScaleMatOperator
  This function multiplies all the current components of the system matrix by the same scalar.

  Collective on Grid

  Input Parameters:
+ grid  - The grid
. alpha - The scalar
- index - The unique operator index returned by GridAddMatOperator()

  Level: intermediate

.keywords grid, matrix, scale
.seealso GridScaleRhs(), GridAddMatOperator()
@*/
int GridScaleMatOperator(Grid grid, PetscScalar alpha, int index)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  grid->matOps[index].alpha *= alpha;
  PetscLogFlops(1);
  PetscFunctionReturn(0);
}