#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: gvec2d.c,v 1.22 2000/10/08 00:27:05 knepley Exp $";
#endif

/* Implements FE vectors derived from 2d triangular grids */
#include "petscsles.h"                 /* For ALE Operators */
#include "src/gvec/gvecimpl.h"         /*I "gvec.h" I*/
#include "src/mesh/impls/triangular/triimpl.h"
#include "gvec2d.h"

#undef  __FUNCT__
#define __FUNCT__ "GVecGetLocalGVec_Triangular_2D"
int GVecGetLocalGVec_Triangular_2D(GVec g, GVec *gvec) {
  SETERRQ(PETSC_ERR_SUP, " ");
}

#undef  __FUNCT__
#define __FUNCT__ "GVecRestoreLocalGVec_Triangular_2D"
int GVecRestoreLocalGVec_Triangular_2D(GVec g, GVec *gvec) {
  SETERRQ(PETSC_ERR_SUP, " ");
}

#undef  __FUNCT__
#define __FUNCT__ "GVecGlobalToLocal_Triangular_2D"
int GVecGlobalToLocal_Triangular_2D(GVec g, InsertMode mode, GVec l) {
  SETERRQ(PETSC_ERR_SUP, " ");
}

#undef  __FUNCT__
#define __FUNCT__ "GVecLocalToGlobal_Triangular_2D"
int GVecLocalToGlobal_Triangular_2D(GVec l, InsertMode mode, GVec g) {
  SETERRQ(PETSC_ERR_SUP, " ");
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateFunction_Triangular_2D"
int GVecEvaluateFunction_Triangular_2D(Grid grid, GVec v, VarOrdering order, PointFunction f, PetscScalar alpha, void *ctx) {
  Mesh          mesh;
  Partition     part;
  int         **localStart = order->localStart;
  FieldClassMap map;
  int           numNodes, numFields;
  int          *fields, **fieldClasses, *classes, *classSizes;
  int           nodeVars;
  PetscScalar  *array;
  double        x, y, z;
  int           size, locSize, overlapSize;
  int           fieldIdx, field, node, nclass, count;
  int           ierr;

  PetscFunctionBegin;
  ierr = GridGetMesh(grid, &mesh);                                                                        CHKERRQ(ierr);
  ierr = MeshGetPartition(mesh, &part);                                                                   CHKERRQ(ierr);
  ierr = VarOrderingGetClassMap(order, &map);                                                             CHKERRQ(ierr);
  numNodes     = map->numNodes;
  numFields    = map->numFields;
  fields       = map->fields;
  fieldClasses = map->fieldClasses;
  classes      = map->classes;
  classSizes   = map->classSizes;
  /* Check for a locally ghosted vector */
  ierr = VecGetArray(v, &array);                                                                          CHKERRQ(ierr);
  ierr = VecGetLocalSize(v, &locSize);                                                                    CHKERRQ(ierr);
  /* ierr = VecGetGhostSize(v, &overlapSize);                                                               CHKERRQ(ierr); */
  overlapSize = locSize + ((Vec_MPI *) v->data)->nghost;
  size        = order->numLocVars;
#if 0
  /* This doesn't work for constrained vector since it conflicts with the check on count
     as variables generated by constraints are not handled here
  */
  if (locSize     != order->numLocVars)
    SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong vector size %d should be %d", locSize, order->numLocVars);
#endif
  if (overlapSize > locSize) {
    ierr = PartitionGetNumOverlapNodes(part, &numNodes);                                                  CHKERRQ(ierr);
    size = order->numOverlapVars;
    if (overlapSize != order->numOverlapVars) SETERRQ(PETSC_ERR_ARG_WRONG, "Wrong size for vector");
  }
  for(node = 0, count = 0; node < numNodes; node++, count += nodeVars) {
    nclass   = classes[node];
    nodeVars = classSizes[nclass];
    for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
      field = fields[fieldIdx];
      if (fieldClasses[fieldIdx][nclass] == 0) continue;
      ierr = MeshGetNodeCoords(mesh, node, &x, &y, &z);                                                   CHKERRQ(ierr);
      ierr = (*f)(1, grid->fields[field].numComp, &x, &y, &z, &array[count+localStart[field][nclass]], ctx); CHKERRQ(ierr);
      array[count+localStart[field][nclass]] *= alpha;
    }
  }
  if (count != size) SETERRQ(PETSC_ERR_PLIB, "Invalid variable offset records");
  ierr = VecRestoreArray(v, &array);                                                                      CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateFunctionBoundary_Triangular_2D"
int GVecEvaluateFunctionBoundary_Triangular_2D(Grid grid, GVec v, int bd, VarOrdering order, PointFunction func,
                                               PetscScalar alpha, void *ctx)
{
  Mesh          mesh         = grid->mesh;
  int         **localStart   = order->localStart;
  int          *offsets      = order->offsets;
  int          *localOffsets = order->localOffsets;
  int           firstVar     = order->firstVar[mesh->part->rank];
  FieldClassMap map;
  int           numNodes, numFields;
  int          *fields, **fieldClasses;
  PetscScalar  *array;        /* The local vector values */
  double        x, y, z;
  int           f, field, node, nclass, row;
  int           ierr;

  PetscFunctionBegin;
  ierr = VarOrderingGetClassMap(order, &map);                                                             CHKERRQ(ierr);
  numNodes     = map->numNodes;
  numFields    = map->numFields;
  fields       = map->fields;
  fieldClasses = map->fieldClasses;
  /* Loop over boundary nodes */
  ierr = VecGetArray(v, &array);                                                                          CHKERRQ(ierr);
  for(f = 0; f < numFields; f++) {
    field = fields[f];
    ierr = (*grid->ops->getboundarystart)(grid, bd, f, PETSC_FALSE, map, &node, &nclass);                 CHKERRQ(ierr);
    while(node >= 0) {
      if (node >= numNodes) {
        row = localOffsets[node-numNodes];
      } else {
        row = offsets[node] - firstVar + localStart[field][nclass];
      }
      ierr = MeshGetNodeCoords(mesh, node, &x, &y, &z);                                                   CHKERRQ(ierr);
      if (fieldClasses[f][nclass] != 0) {
        ierr = (*func)(1, grid->fields[field].numComp, &x, &y, &z, &array[row], ctx);                     CHKERRQ(ierr);
        array[row] *= alpha;
      }
      ierr = (*grid->ops->getboundarynext)(grid, bd, f, PETSC_FALSE, map, &node, &nclass);                CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                       CHKERRQ(ierr);
#endif
    }
  }
  ierr = VecRestoreArray(v, &array);                                                                      CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateFunctionCollective_Triangular_2D"
int GVecEvaluateFunctionCollective_Triangular_2D(Grid grid, GVec v, VarOrdering order, PointFunction f, PetscScalar alpha,
                                                 void *ctx)
{
  Mesh             mesh;
  FieldClassMap    map;
  int            **localStart = order->localStart;
  int              numNodes, numFields;
  int             *fields, **fieldClasses, *classes, *classSizes;
  int              nodeVars, comp;
  PetscScalar     *array;
  double           x, y, z;
  int              maxNodes; /* The most nodes in any domain */
  int              fieldIdx, field, node, nclass, count;
  int              ierr;

  PetscFunctionBegin;
  ierr = GridGetMesh(grid, &mesh);                                                                        CHKERRQ(ierr);
  ierr = VarOrderingGetClassMap(order, &map);                                                             CHKERRQ(ierr);
  numNodes     = map->numNodes;
  numFields    = map->numFields;
  fields       = map->fields;
  fieldClasses = map->fieldClasses;
  classes      = map->classes;
  classSizes   = map->classSizes;
  ierr = MPI_Allreduce(&numNodes, &maxNodes, 1, MPI_INT, MPI_MAX, grid->comm);                            CHKERRQ(ierr);
  ierr = VecGetArray(v, &array);                                                                          CHKERRQ(ierr);
  for(node = 0, count = 0; node < maxNodes; node++) {
    if (node < numNodes) {
      nclass     = classes[node];
      nodeVars   = classSizes[nclass];
      for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
        field = fields[fieldIdx];
        comp  = grid->fields[field].numComp;
        if (fieldClasses[fieldIdx][nclass] == 0) {
          /* We have to make sure that every processor is available at each iteration */
          ierr = (*f)(0, 0, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, ctx);                         CHKERRQ(ierr);
          continue;
        }
        ierr = MeshGetNodeCoords(mesh, node, &x, &y, &z);                                                 CHKERRQ(ierr);
        ierr = (*f)(1, comp, &x, &y, &z, &array[count+localStart[field][nclass]], ctx);                   CHKERRQ(ierr);
        array[count+localStart[field][nclass]] *= alpha;
      }
      count += nodeVars;
    } else {
      /* We have to make sure that every processor is available at each iteration */
      ierr = (*f)(0, 0, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL, ctx);                             CHKERRQ(ierr);
    }
  }
  if (count != order->numLocVars) {
    SETERRQ2(PETSC_ERR_PLIB, "Invalid number of variables modified %d should be %d", count, order->numLocVars);
  }
  ierr = VecRestoreArray(v, &array);                                                                      CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateFunctionGalerkin_Triangular_2D"
int GVecEvaluateFunctionGalerkin_Triangular_2D(Grid grid, GVec v, int numFields, int *fields, LocalVarOrdering locOrder,
                                               PointFunction f, PetscScalar alpha, void *ctx)
{
  Mesh         mesh        = grid->mesh;
  int          numElements = mesh->numFaces;
  ElementVec   vec         = grid->vec;
  int         *elemStart   = locOrder->elemStart;
  PetscScalar *array       = vec->array;
  int          field, fieldIdx, elem;
#ifdef PETSC_USE_BOPT_g
  PetscTruth   opt;
#endif
  int          ierr;

  PetscFunctionBegin;
  /* Loop over elements */
  for(elem = 0; elem < numElements; elem++) {
    /* Initialize element vector */
    ierr = ElementVecZero(vec);                                                                           CHKERRQ(ierr);

    /* Get contribution to the element vector from each discretization */
    for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
      field = fields[fieldIdx];
      ierr  = DiscretizationEvaluateFunctionGalerkin(grid->fields[field].disc, mesh, f, alpha, elem, &array[elemStart[field]], ctx);
      CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                        CHKERRQ(ierr);
#endif
    }

    /* Setup global row and column indices */
    ierr = GridCalcElementVecIndices(grid, elem, vec);                                                    CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
    ierr = PetscOptionsHasName(PETSC_NULL, "-trace_vec_assembly", &opt);                                  CHKERRQ(ierr);
    if (opt == PETSC_TRUE) {
      int var;

      for(var = 0; var < vec->reduceSize; var++)
        PetscPrintf(PETSC_COMM_SELF, "%2d %4.2g\n", vec->indices[var], PetscRealPart(array[var]));
    }
#endif
    /* Put values in global vector */
    ierr = ElementVecSetValues(vec, v, ADD_VALUES);                                                       CHKERRQ(ierr);
  }

  ierr = VecAssemblyBegin(v);                                                                             CHKERRQ(ierr);
  ierr = VecAssemblyEnd(v);                                                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateFunctionGalerkinCollective_Triangular_2D"
int GVecEvaluateFunctionGalerkinCollective_Triangular_2D(Grid grid, GVec v, int numFields, int *fields,
                                                         LocalVarOrdering locOrder, PointFunction f,
                                                         PetscScalar alpha, void *ctx)
{
  Mesh         mesh        = grid->mesh;
  int          numElements = mesh->numFaces;
  ElementVec   vec         = grid->vec;
  int         *elemStart   = locOrder->elemStart;
  PetscScalar *array       = vec->array;
  int          maxElements;
  int          field, fieldIdx, elem;
#ifdef PETSC_USE_BOPT_g
  PetscTruth   opt;
#endif
  int          ierr;

  PetscFunctionBegin;
  ierr = MPI_Allreduce(&numElements, &maxElements, 1, MPI_INT, MPI_MAX, grid->comm);                     CHKERRQ(ierr);
  /* Loop over elements */
  for(elem = 0; elem < maxElements; elem++) {
    if (elem < numElements) {
      /* Initialize element vector */
      ierr = ElementVecZero(vec);                                                                         CHKERRQ(ierr);

      /* Get contribution to the element vector from each discretization */
      for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
        field = fields[fieldIdx];
        ierr  = DiscretizationEvaluateFunctionGalerkin(grid->fields[field].disc, mesh, f, alpha, elem, &array[elemStart[field]], ctx);
        CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
        ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                      CHKERRQ(ierr);
#endif
      }

      /* Setup global row and column indices */
      ierr = GridCalcElementVecIndices(grid, elem, vec);                                                  CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscOptionsHasName(PETSC_NULL, "-trace_vec_assembly", &opt);                                CHKERRQ(ierr);
      if (opt == PETSC_TRUE) {
        int var;

        for(var = 0; var < vec->reduceSize; var++)
          PetscPrintf(PETSC_COMM_SELF, "%2d %4.2g\n", vec->indices[var], PetscRealPart(array[var]));
      }
#endif
      /* Put values in global vector */
      ierr = ElementVecSetValues(vec, v, ADD_VALUES);                                                     CHKERRQ(ierr);
    } else {
      /* We have to make sure that every processor is available at each call to f */
      for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
        field = fields[fieldIdx];
        ierr  = DiscretizationEvaluateFunctionGalerkin(grid->fields[field].disc, mesh, f, 0.0, -1, PETSC_NULL, ctx);   CHKERRQ(ierr);
      }
    }
  }

  ierr = VecAssemblyBegin(v);                                                                             CHKERRQ(ierr);
  ierr = VecAssemblyEnd(v);                                                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateBoundaryFunctionGalerkin_Triangular_2D"
int GVecEvaluateBoundaryFunctionGalerkin_Triangular_2D(Grid grid, GVec v, int numFields, int *fields,
                                                       LocalVarOrdering locOrder, PointFunction f, PetscScalar alpha, void *ctx)
{
  Mesh                     mesh      = grid->mesh;
  Partition                part;
  Mesh_Triangular         *tri       = (Mesh_Triangular *) grid->mesh->data;
  int                      elemSize  = locOrder->elemSize;
  int                     *elemStart = locOrder->elemStart;
  int                      numEdges  = mesh->numEdges;
  int                     *bdEdges   = tri->bdEdges;
  int                      firstEdge;
  ElementVec               vec;         /* The element vector */
  PetscScalar             *array;       /* The values in the element vector */
  EdgeContext              bdCtx;       /* A context wrapper to communicate the midnode of an edge */
  int                      field, edge, midNode;
  int                      fieldIdx, bd, bdEdge;
#ifdef PETSC_USE_BOPT_g
  PetscTruth               opt;
#endif
  int                      ierr;

  PetscFunctionBegin;
  /* Setup element vector for the lower dimensional system */
  ierr  = ElementVecCreate(grid->comm, elemSize, &vec);                                                   CHKERRQ(ierr);
  array = vec->array;

  /* Setup user context */
  bdCtx.ctx = ctx;

  /* Our problem here is that "edges" are not data structures like "elements". The element
     holds the midnodes which appear on it, but edges do not. Thus we must pass the midnode
     number to the discretization, which we do using a context wrapper. Unfortunately, the
     row indices were derived from element, so we must introduce another numbering function
     which operates on nodes alone. The midnode number is found by a search of the elements
     which could certainly be improved with geometric hints. We might also assume that it
     is the node lying between the two endpoints in the bdNodes[] array. In addition, the
     boundary variable ordering is in relation to boundary node numbers, so that the node
     number must be converted before calling the numbering function. This could be speeded up
     by placing boundary node numbers in the bdEdges[] array instead. */

  /* Loop over boundary edges */
  ierr = MeshGetPartition(mesh, &part);                                                                   CHKERRQ(ierr);
  ierr = PartitionGetStartEdge(part, &firstEdge);                                                         CHKERRQ(ierr);
  for(bd = 0, bdEdge = 0; bd < grid->numBd; bd++) {
    for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) {
      /* Check that edge is on this processor */
      edge = bdEdges[bdEdge] - firstEdge;
      if ((edge < 0) || (edge > numEdges)) continue;

      /* Search for midnode on edge */
      midNode = -1;
      ierr = MeshGetMidnodeFromEdge(mesh, edge, &midNode);                                                CHKERRQ(ierr);
      bdCtx.midnode = midNode;

      /* Initialize element matrix */
      ierr = ElementVecZero(vec);                                                                         CHKERRQ(ierr);

      for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
        field = fields[fieldIdx];
        ierr = DiscretizationEvaluateFunctionGalerkin(grid->fields[field].disc->bdDisc, mesh, f, alpha, edge, &array[elemStart[field]], &bdCtx);
        CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
        ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                      CHKERRQ(ierr);
#endif
      }

      /* Setup global row and column indices */
      ierr = GridCalcBoundaryElementVecIndices(grid, bd, edge, midNode, grid->bdOrder, PETSC_FALSE, vec); CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscOptionsHasName(PETSC_NULL, "-trace_vec_bd_assembly", &opt);                             CHKERRQ(ierr);
      if (opt == PETSC_TRUE) {
        int var;

        for(var = 0; var < vec->reduceSize; var++)
          PetscPrintf(PETSC_COMM_SELF, "%2d %4.2g\n", vec->indices[var], PetscRealPart(array[var]));
      }
#endif
      /* Put values in global vector */
      ierr = ElementVecSetValues(vec, v, ADD_VALUES);                                                     CHKERRQ(ierr);
    }
  }
#ifdef PETSC_USE_BOPT_g
  if (bdEdge != mesh->numBdEdges) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary edge numbering");
#endif

  ierr = VecAssemblyBegin(v);                                                                             CHKERRQ(ierr);
  ierr = VecAssemblyEnd(v);                                                                               CHKERRQ(ierr);

  /* Cleanup */
  ierr = ElementVecDestroy(vec);                                                                          CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateBoundaryFunctionGalerkinCollective_Triangular_2D"
int GVecEvaluateBoundaryFunctionGalerkinCollective_Triangular_2D(Grid grid, GVec v, int numFields, int *fields,
                                                                 LocalVarOrdering locOrder, PointFunction f,
                                                                 PetscScalar alpha, void *ctx)
{
  Mesh                     mesh      = grid->mesh;
  Partition                part;
  Mesh_Triangular         *tri       = (Mesh_Triangular *) mesh->data;
  int                      elemSize  = locOrder->elemSize;
  int                     *elemStart = locOrder->elemStart;
  int                      numEdges  = mesh->numEdges;
  int                     *bdEdges   = tri->bdEdges;
  int                      firstEdge;
  ElementVec               vec;         /* The element vector */
  PetscScalar             *array;       /* The values in the element vector */
  EdgeContext              bdCtx;       /* A context wrapper to communicate the midnode of an edge */
  int                      field, edge, midnode;
  int                      fieldIdx, bd, bdEdge;
#ifdef PETSC_USE_BOPT_g
  PetscTruth               opt;
#endif
  int                      ierr;

  PetscFunctionBegin;
  /* Setup element vector for the lower dimensional system */
  ierr  = ElementVecCreate(grid->comm, elemSize, &vec);                                                   CHKERRQ(ierr);
  array = vec->array;

  /* Setup user context */
  bdCtx.ctx = ctx;

  /* Our problem here is that "edges" are not data structures like "elements". The element
     holds the midnodes which appear on it, but edges do not. Thus we must pass the midnode
     number to the discretization, which we do using a context wrapper. Unfortunately, the
     row indices were derived from element, so we must introduce another numbering function
     which operates on nodes alone. The midnode number is found by a search of the elements
     which could certainly be improved with geometric hints. We might also assume that it
     is the node lying between the two endpoints in the bdNodes[] array. In addition, the
     boundary variable ordering is in relation to boundary node numbers, so that the node
     number must be converted before calling the numbering function. This could be speeded up
     by placing boundary node numbers in the bdEdges[] array instead. */

  /* Loop over boundary edges */
  ierr = MeshGetPartition(mesh, &part);                                                                   CHKERRQ(ierr);
  ierr = PartitionGetStartEdge(part, &firstEdge);                                                         CHKERRQ(ierr);
  for(bd = 0, bdEdge = 0; bd < grid->numBd; bd++) {
    for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) {
      /* Check that edge is on this processor */
      edge = bdEdges[bdEdge] - firstEdge;
      if ((edge < 0) || (edge > numEdges)) {
        for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
          field = fields[fieldIdx];
          ierr  = DiscretizationEvaluateFunctionGalerkin(grid->fields[field].disc->bdDisc, mesh, f, 0.0, -1, PETSC_NULL, &bdCtx);
          CHKERRQ(ierr);
        }
        continue;
      }

      /* Locate midnode on edge */
      midnode = -1;
      ierr = MeshGetMidnodeFromEdge(mesh, edge, &midnode);                                                CHKERRQ(ierr);
      bdCtx.midnode = midnode;
#ifdef PETSC_USE_BOPT_g
      if (tri->markers[midnode] != tri->bdMarkers[bd])
        SETERRQ4(PETSC_ERR_ARG_WRONG, "Invalid midnode %d has marker %d on boundary %d (%d)",
                 midnode, tri->markers[midnode], bd, tri->bdMarkers[bd]);
#endif

      /* Initialize element matrix */
      ierr = ElementVecZero(vec);                                                                         CHKERRQ(ierr);

      for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
        field = fields[fieldIdx];
        ierr = DiscretizationEvaluateFunctionGalerkin(grid->fields[field].disc->bdDisc, mesh, f, alpha, edge, &array[elemStart[field]], &bdCtx);
        CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
        ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                      CHKERRQ(ierr);
#endif
      }

      /* Setup global row and column indices */
      ierr = GridCalcBoundaryElementVecIndices(grid, bd, edge, midnode, grid->bdOrder, PETSC_FALSE, vec); CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscOptionsHasName(PETSC_NULL, "-trace_vec_bd_assembly", &opt);                             CHKERRQ(ierr);
      if (opt == PETSC_TRUE) {
        int var;

        for(var = 0; var < vec->reduceSize; var++)
          PetscPrintf(PETSC_COMM_SELF, "%2d %4.2g\n", vec->indices[var], PetscRealPart(array[var]));
      }
#endif
      /* Put values in global vector */
      ierr = ElementVecSetValues(vec, v, ADD_VALUES);                                                     CHKERRQ(ierr);
    }
  }
#ifdef PETSC_USE_BOPT_g
  if (bdEdge != mesh->numBdEdges) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary edge numbering");
#endif

  ierr = VecAssemblyBegin(v);                                                                             CHKERRQ(ierr);
  ierr = VecAssemblyEnd(v);                                                                               CHKERRQ(ierr);

  /* Cleanup */
  ierr = ElementVecDestroy(vec);                                                                          CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateNonlinearOperatorGalerkin_Triangular_2D"
int GVecEvaluateNonlinearOperatorGalerkin_Triangular_2D(Grid grid, GVec v, GVec x, GVec y, int numFields, int *fields,
                                                        LocalVarOrdering locOrder, NonlinearOperator op, PetscScalar alpha,
                                                        PetscTruth isALE, void *ctx)
{
  Mesh         mesh;
  Partition    part;
  MeshMover    mover;
  int         *elemStart     = locOrder->elemStart;
  ElementVec   vec           = grid->vec;
  PetscScalar *array         = vec->array;
  ElementVec   ghostVec      = grid->ghostElementVec; /* The local solution vector */
  PetscScalar *ghostArray    = ghostVec->array;       /* The values in the ghost element vector */
  PetscTruth   reduceElement = grid->reduceElementArgs;
  Grid         ALEGrid;                             /* The grid describing the mesh velocity */
  Vec          appVec;                              /* The local vec for y */
  ElementVec   elemAppVec;                          /* The element vec for y */
  PetscScalar *appArray;                            /* The values in elemAppVec */
  ElementVec   MeshALEVec;                          /* ALE velocity vector with mesh discretization */
  ElementVec   ALEVec;                              /* ALE velocity vector */
  PetscScalar *ALEArray;                            /* The values in the ALE element vector */
  PetscScalar *nonlinearArgs[2];
  PetscTruth   ALEActive;
  int          numElements;
  int          field, fieldIdx, elem;
#ifdef PETSC_USE_BOPT_g
  PetscTruth   opt;
#endif
  int          ierr;

  PetscFunctionBegin;
  ierr = GridGetMesh(grid, &mesh);                                                                        CHKERRQ(ierr);
  ierr = MeshGetPartition(mesh, &part);                                                                   CHKERRQ(ierr);
  if (grid->ALEActive && (isALE == PETSC_TRUE)) {
    ALEActive = PETSC_TRUE;
    ierr = MeshGetMover(mesh, &mover);                                                                    CHKERRQ(ierr);
    ierr = MeshMoverGetVelocityGrid(mover, &ALEGrid);                                                     CHKERRQ(ierr);
  } else {
    ALEActive = PETSC_FALSE;
  }
  /* Fill the local solution vectors */
  if (x != PETSC_NULL) {
    ierr = GridGlobalToLocal(grid, INSERT_VALUES, x);                                                     CHKERRQ(ierr);
  }
  ierr = VecDuplicate(grid->ghostVec, &appVec);                                                           CHKERRQ(ierr);
  ierr = ElementVecDuplicate(ghostVec, &elemAppVec);                                                      CHKERRQ(ierr);
  if (y != PETSC_NULL) {
    ierr = GridGlobalToLocalGeneral(grid, y, appVec, INSERT_VALUES, grid->ghostScatter);                  CHKERRQ(ierr);
  }
  appArray = elemAppVec->array;

  /* Setup ALE variables */
  if (ALEActive == PETSC_TRUE) {
    /* Notice that the ALEArray is from this grid, not the mesh velocity grid */
    MeshALEVec = ALEGrid->vec;
    ierr       = ElementVecDuplicate(grid->vec, &ALEVec);                                                 CHKERRQ(ierr);
    ALEArray   = ALEVec->array;
  } else {
    MeshALEVec = PETSC_NULL;
    ALEArray   = PETSC_NULL;
  }

  /* Loop over elements */
  ierr = PartitionGetNumElements(part, &numElements);                                                     CHKERRQ(ierr);
  for(elem = 0; elem < numElements; elem++) {
    /* Initialize element vector */
    ierr = ElementVecZero(vec);                                                                           CHKERRQ(ierr);

    /* Setup local row and column indices */
    ierr = GridCalcLocalElementVecIndices(grid, elem, ghostVec);                                          CHKERRQ(ierr);
    ierr = ElementVecDuplicateIndices(ghostVec, elemAppVec);                                              CHKERRQ(ierr);

    /* Setup local solution vector */
    ierr = GridLocalToElement(grid, ghostVec);                                                            CHKERRQ(ierr);
    ierr = GridLocalToElementGeneral(grid, appVec, grid->bdReduceVecCur, grid->reduceSystem, reduceElement, elemAppVec);CHKERRQ(ierr);

    /* Setup ALE variables */
    if (ALEActive == PETSC_TRUE) {
      ierr = GridCalcLocalElementVecIndices(ALEGrid, elem, MeshALEVec);                                   CHKERRQ(ierr);
      ierr = GridLocalToElement(ALEGrid, MeshALEVec);                                                     CHKERRQ(ierr);
    }

    /* Get contribution to the element vector from each discretization */
    for(fieldIdx = 0; fieldIdx < numFields; fieldIdx++) {
      field            = fields[fieldIdx];
      nonlinearArgs[0] = &ghostArray[elemStart[field]];
      nonlinearArgs[1] = &appArray[elemStart[field]];
      if (ALEActive == PETSC_TRUE)
      {
        ierr = GridInterpolateElementVec(ALEGrid, 0, MeshALEVec, grid, field, ALEVec);                    CHKERRQ(ierr);
        ierr = DiscretizationEvaluateNonlinearALEOperatorGalerkin(grid->fields[field].disc, mesh, op, alpha, elem, 2, nonlinearArgs,
                                                                  ALEArray, &array[elemStart[field]], ctx);
        CHKERRQ(ierr);
      } else {
        ierr = DiscretizationEvaluateNonlinearOperatorGalerkin(grid->fields[field].disc, mesh, op, alpha, elem, 2, nonlinearArgs,
                                                               &array[elemStart[field]], ctx);
        CHKERRQ(ierr);
      }
#ifdef PETSC_USE_BOPT_g
      ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                        CHKERRQ(ierr);
#endif
    }

    /* Setup global row and column indices */
    ierr = GridCalcElementVecIndices(grid, elem, vec);                                                    CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
    ierr = PetscOptionsHasName(PETSC_NULL, "-trace_vec_assembly", &opt);                                  CHKERRQ(ierr);
    if (opt == PETSC_TRUE) {
      int var;

      for(var = 0; var < vec->reduceSize; var++)
        PetscPrintf(PETSC_COMM_SELF, "%2d %4.2g\n", vec->indices[var], PetscRealPart(array[var]));
    }
#endif
    /* Put values in global vector */
    ierr = ElementVecSetValues(vec, v, ADD_VALUES);                                                       CHKERRQ(ierr);
  }

  /* Cleanup ALE variables */
  if (ALEActive == PETSC_TRUE) {
    ierr = ElementVecDestroy(ALEVec);                                                                     CHKERRQ(ierr);
  }
  /* Cleanup additional input vectors */
  ierr = VecDestroy(appVec);                                                                              CHKERRQ(ierr);
  ierr = ElementVecDestroy(elemAppVec);                                                                   CHKERRQ(ierr);
  ierr = VecAssemblyBegin(v);                                                                             CHKERRQ(ierr);
  ierr = VecAssemblyEnd(v);                                                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateOperatorGalerkin_Triangular_2D"
int GVecEvaluateOperatorGalerkin_Triangular_2D(Grid grid, GVec v, GVec x, GVec y, VarOrdering sOrder, LocalVarOrdering sLocOrder,
                                               VarOrdering tOrder, LocalVarOrdering tLocOrder, int op, PetscScalar alpha, void *ctx)
{
  Mesh             mesh          = grid->mesh;
  PetscTruth       reduceSystem  = grid->reduceSystem;
  PetscTruth       reduceElement = grid->reduceElement;
  int              sElemSize     = sLocOrder->elemSize;
  int             *sElemStart    = sLocOrder->elemStart;
  int              tElemSize     = tLocOrder->elemSize;
  int             *tElemStart    = tLocOrder->elemStart;
  FieldClassMap    sMap,         tMap;
  int              numSFields,   numTFields;
  int             *sFields,     *tFields;
  PetscTruth       sConstrained, tConstrained;
  Vec              ghostVec;     /* The local ghost vector for x (usually the solution) */
  VecScatter       ghostScatter; /* The scatter from x to ghostVec */
  Vec              appVec;       /* The local ghost vector for y (usually the application vector) */
  VecScatter       appScatter;   /* The scatter from y to appVec */
  ElementMat       mat;
  ElementVec       elemGhostVec, elemAppVec, vec;
  PetscScalar     *ghostArray, *appArray, *matArray, *array;
  int              numElements;
  int              f, sField, tField, elem;
#ifdef PETSC_USE_BOPT_g
  PetscTruth            opt;
#endif
  int              ierr;

  PetscFunctionBegin;
  ierr = MeshGetInfo(mesh, PETSC_NULL, PETSC_NULL, PETSC_NULL, &numElements);                             CHKERRQ(ierr);
  ierr = VarOrderingGetClassMap(sOrder, &sMap);                                                           CHKERRQ(ierr);
  ierr = VarOrderingGetClassMap(tOrder, &tMap);                                                           CHKERRQ(ierr);
  numSFields   = sMap->numFields;
  sFields      = sMap->fields;
  sConstrained = sMap->isConstrained;
  numTFields   = tMap->numFields;
  tFields      = tMap->fields;
  tConstrained = tMap->isConstrained;
  /* Setup reduction */
  ierr = (*grid->ops->gridsetupghostscatter)(grid, tOrder, &ghostVec, &ghostScatter);                     CHKERRQ(ierr);
  ierr = (*grid->ops->gridsetupghostscatter)(grid, sOrder, &appVec,   &appScatter);                       CHKERRQ(ierr);
  /* Setup element vector and matrix */
  if (tConstrained == PETSC_TRUE) {
    for(f = 0; f < numTFields; f++) {
      if (grid->fields[tFields[f]].isConstrained == PETSC_TRUE)
        tElemSize += grid->fields[tFields[f]].disc->funcs*grid->fields[tFields[f]].constraintCompDiff;
    }
  }
  if (sConstrained == PETSC_TRUE) {
    for(f = 0; f < numSFields; f++) {
      if (grid->fields[sFields[f]].isConstrained == PETSC_TRUE)
        sElemSize += grid->fields[sFields[f]].disc->funcs*grid->fields[sFields[f]].constraintCompDiff;
    }
  }
  ierr       = ElementVecCreate(grid->comm, tElemSize, &vec);                                             CHKERRQ(ierr);
  array      = vec->array;
  ierr       = ElementVecDuplicate(vec, &elemGhostVec);                                                   CHKERRQ(ierr);
  ghostArray = elemGhostVec->array;
  ierr       = ElementVecCreate(grid->comm, sElemSize, &elemAppVec);                                      CHKERRQ(ierr);
  appArray   = elemAppVec->array;
  ierr       = ElementMatCreate(grid->comm, tElemSize, sElemSize, &mat);                                  CHKERRQ(ierr);
  matArray   = mat->array;

  /* Fill the local solution vectors */
  ierr = GridGlobalToLocalGeneral(grid, x, ghostVec, INSERT_VALUES, ghostScatter);                        CHKERRQ(ierr);
  ierr = GridGlobalToLocalGeneral(grid, y, appVec,   INSERT_VALUES, appScatter);                          CHKERRQ(ierr);

  /* Setup the operator with information about the test function space */
  for(f = 0; f < numSFields; f++) {
    grid->fields[sFields[f]].disc->operators[op]->test = grid->fields[tFields[f]].disc;
  }

  /* Loop over elements */
  for(elem = 0; elem < numElements; elem++) {
    /* Initialize element vector */
    ierr = ElementVecZero(vec);                                                                           CHKERRQ(ierr);
    vec->reduceSize          = tLocOrder->elemSize;
    elemGhostVec->reduceSize = tLocOrder->elemSize;
    elemAppVec->reduceSize   = sLocOrder->elemSize;

    /* Setup local row indices */
    ierr = GridCalcGeneralElementVecIndices(grid, elem, tOrder, PETSC_NULL, PETSC_TRUE, elemGhostVec);    CHKERRQ(ierr);
    ierr = GridCalcGeneralElementVecIndices(grid, elem, sOrder, PETSC_NULL, PETSC_TRUE, elemAppVec);      CHKERRQ(ierr);
    /* Setup local vectors */
    ierr = GridLocalToElementGeneral(grid, ghostVec, grid->bdReduceVecCur, reduceSystem, reduceElement, elemGhostVec);CHKERRQ(ierr);
    ierr = GridLocalToElementGeneral(grid, appVec,   grid->bdReduceVecCur, reduceSystem, reduceElement, elemAppVec);CHKERRQ(ierr);
    /* Must transform to unconstrained variables for element integrals */
    ierr = GridProjectElementVec(grid, mesh, elem, tOrder, PETSC_FALSE, elemGhostVec);                    CHKERRQ(ierr);
    ierr = GridProjectElementVec(grid, mesh, elem, sOrder, PETSC_FALSE, elemAppVec);                      CHKERRQ(ierr);

    for(f = 0; f < numSFields; f++) {
      sField = sFields[f];
      tField = tFields[f];
      /* Get contribution to the element vector from the linear operator */
      ierr = ElementMatZero(mat);                                                                         CHKERRQ(ierr);
      ierr = DiscretizationEvaluateOperatorGalerkinMF(grid->fields[sField].disc, mesh, sElemSize, tElemStart[tField], sElemStart[sField],
                                                      op, alpha, elem, &ghostArray[sElemStart[sField]],
                                                      &appArray[sElemStart[sField]], array, matArray, ctx);
      CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                       CHKERRQ(ierr);
#endif
    }

    /* Setup global row indices, with reduction if necessary */
    ierr = GridCalcGeneralElementVecIndices(grid, elem, tOrder, PETSC_NULL, PETSC_FALSE, vec);            CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
    ierr = PetscOptionsHasName(PETSC_NULL, "-trace_vec_assembly", &opt);                                  CHKERRQ(ierr);
    if (opt == PETSC_TRUE) {
      int var;

      for(var = 0; var < vec->reduceSize; var++)
        PetscPrintf(PETSC_COMM_SELF, "%2d %4.2g\n", vec->indices[var], PetscRealPart(array[var]));
    }
#endif
    /* Put values in global vector */
    ierr = ElementVecSetValues(vec, v, ADD_VALUES);                                                       CHKERRQ(ierr);
  }

  ierr = VecDestroy(ghostVec);                                                                            CHKERRQ(ierr);
  ierr = VecScatterDestroy(ghostScatter);                                                                 CHKERRQ(ierr);
  ierr = VecDestroy(appVec);                                                                              CHKERRQ(ierr);
  ierr = VecScatterDestroy(appScatter);                                                                   CHKERRQ(ierr);
  ierr = ElementVecDestroy(elemGhostVec);                                                                 CHKERRQ(ierr);
  ierr = ElementVecDestroy(elemAppVec);                                                                   CHKERRQ(ierr);
  ierr = ElementVecDestroy(vec);                                                                          CHKERRQ(ierr);
  ierr = ElementMatDestroy(mat);                                                                          CHKERRQ(ierr);
  ierr = VecAssemblyBegin(v);                                                                             CHKERRQ(ierr);
  ierr = VecAssemblyEnd(v);                                                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateSystemMatrix_Triangular_2D"
int GVecEvaluateSystemMatrix_Triangular_2D(Grid grid, GVec x, GVec y, GVec f, void *ctx)
{
  Mesh                  mesh                = grid->mesh;
  int                   numElements         = mesh->numFaces;
  int                   numMatOps           = grid->numMatOps;       /* The number of operators in the matrix */
  GridOp               *matOps              = grid->matOps;          /* The operators in the system matrix */
  PetscTruth            reduceSystem        = grid->reduceSystem;
  PetscTruth            reduceElement       = grid->reduceElement;
  PetscTruth            explicitConstraints = grid->explicitConstraints;
  PetscConstraintObject constCtx            = grid->constraintCtx;   /* The constraint object */
  int                   numNewFields        = grid->numNewFields;    /* The number of new fields added by constraints */
  ElementVec            vec                 = grid->vec;             /* The element vector */
  PetscScalar          *array               = vec->array;            /* The values in the element vector */
  ElementMat            mat                 = grid->mat;             /* The element matrix */
  PetscScalar          *matArray            = mat->array;            /* The values in the element matrix */
  Vec                   ghostVec            = grid->ghostVec;        /* The local solution vector */
  ElementVec            elemGhostVec        = grid->ghostElementVec; /* Local solution vector */
  PetscScalar          *ghostArray          = elemGhostVec->array;   /* The values in the ghost element vector */
  int                   numFields           = grid->cm->numFields;   /* The number of fields in the calculation */
  int                  *fields              = grid->cm->fields;      /* The fields participating in the calculation */
  LocalVarOrdering      locOrder            = grid->locOrder;        /* The default local variable ordering */
  int                   elemSize            = locOrder->elemSize;    /* The number of shape functions in the element matrix */
  int                  *elemStart           = locOrder->elemStart;   /* The offset of each field in the element matrix */
  int                   rank                = mesh->part->rank;      /* The processor rank */
  MeshMover             mover;
  Grid                  ALEGrid;                                     /* The grid describing the mesh velocity */
  VarOrdering           order;                                       /* The default variable ordering */
  ElementVec            MeshALEVec;                                  /* ALE velocity vector with mesh discretization */
  ElementVec            ALEVec;                                      /* ALE velocity vector */
  PetscScalar          *ALEArray;                                    /* The values in the ALE element vector */
  Vec                   appVec;                                      /* The local vec for y */
  ElementVec            elemAppVec;                                  /* The element vec for y */
  PetscScalar          *appArray;                                    /* The values in elemAppVec */
  int                   elem, fieldIndex;
  int                   newComp = 0;
  int                   sField, tField, op, newField, row, col;
#ifdef PETSC_USE_BOPT_g
  PetscTruth            opt;
#endif
  int                   ierr;

  PetscFunctionBegin;
  ierr = MeshGetMover(mesh, &mover);                                                                      CHKERRQ(ierr);
  ierr = MeshMoverGetVelocityGrid(mover, &ALEGrid);                                                       CHKERRQ(ierr);
  ierr = PetscObjectQuery((PetscObject) x, "Order", (PetscObject *) &order);                              CHKERRQ(ierr);
  /* Right now, we ignore the preconditioner */
  /* Fill the local solution vectors */
  if (x != PETSC_NULL) {
    ierr = GridGlobalToLocal(grid, INSERT_VALUES, x);                                                     CHKERRQ(ierr);
  }
  ierr = VecDuplicate(ghostVec, &appVec);                                                                 CHKERRQ(ierr);
  ierr = GridGlobalToLocalGeneral(grid, y, appVec, INSERT_VALUES, grid->ghostScatter);                    CHKERRQ(ierr);
  ierr = ElementVecDuplicate(elemGhostVec, &elemAppVec);                                                  CHKERRQ(ierr);
  appArray = elemAppVec->array;

  /* Setup ALE variables */
  if (grid->ALEActive == PETSC_TRUE) {
    /* Notice that the ALEArray is from this grid, not the mesh velocity grid */
    MeshALEVec = ALEGrid->vec;
    ierr       = ElementVecDuplicate(grid->vec, &ALEVec);                                                 CHKERRQ(ierr);
    ALEArray   = ALEVec->array;
  } else {
    MeshALEVec = PETSC_NULL;
    ALEArray   = PETSC_NULL;
  }

  /* Loop over elements */
  for(elem = 0; elem < numElements; elem++)
  {
    /* Initialize element vector */
    ierr = ElementVecZero(vec);                                                                           CHKERRQ(ierr);
    vec->reduceSize          = locOrder->elemSize;
    elemGhostVec->reduceSize = locOrder->elemSize;
    elemAppVec->reduceSize   = locOrder->elemSize;

    /* Setup global row and column indices */
    ierr = GridCalcLocalElementVecIndices(grid, elem, elemGhostVec);                                      CHKERRQ(ierr);
    ierr = ElementVecDuplicateIndices(elemGhostVec, elemAppVec);                                          CHKERRQ(ierr);
    elemAppVec->reduceSize = elemGhostVec->reduceSize;

    /* Setup local solution vector */
    ierr = GridLocalToElementGeneral(grid, ghostVec, grid->bdReduceVecCur, reduceSystem, reduceElement, elemGhostVec);CHKERRQ(ierr);
    ierr = GridLocalToElementGeneral(grid, appVec, grid->bdReduceVecCur, reduceSystem, reduceElement, elemAppVec);CHKERRQ(ierr);

    /* Must transform to unconstrained variables for element integrals */
    ierr = GridProjectElementVec(grid, mesh, elem, order, PETSC_FALSE, elemGhostVec);                     CHKERRQ(ierr);
    ierr = GridProjectElementVec(grid, mesh, elem, order, PETSC_FALSE, elemAppVec);                       CHKERRQ(ierr);

    /* Setup ALE variables */
    if (grid->ALEActive == PETSC_TRUE) {
      ierr = GridCalcLocalElementVecIndices(ALEGrid, elem, MeshALEVec);                                   CHKERRQ(ierr);
      ierr = GridLocalToElement(ALEGrid, MeshALEVec);                                                     CHKERRQ(ierr);
    }

    /* Calculate the contribution to the element matrix from each field */
    for(op = 0; op < numMatOps; op++) {
      sField = matOps[op].field;
      tField = grid->fields[sField].disc->operators[matOps[op].op]->test->field;
      if (grid->fields[sField].isActive) {
        ierr = ElementMatZero(mat);                                                                       CHKERRQ(ierr);
        if (matOps[op].isALE) {
          ierr = GridInterpolateElementVec(ALEGrid, 0, MeshALEVec, grid, sField, ALEVec);                 CHKERRQ(ierr);
          ierr = DiscretizationEvaluateALEOperatorGalerkinMF(grid->fields[sField].disc, mesh, elemSize, elemStart[tField], elemStart[sField],
                                                             matOps[op].op, matOps[op].alpha, elem, &ghostArray[elemStart[sField]],
                                                             &appArray[elemStart[sField]], ALEArray, array, matArray, ctx);
          CHKERRQ(ierr);
        } else {
          ierr = DiscretizationEvaluateOperatorGalerkinMF(grid->fields[sField].disc, mesh, elemSize, elemStart[tField], elemStart[sField],
                                                          matOps[op].op, matOps[op].alpha, elem, &ghostArray[elemStart[sField]],
                                                          &appArray[elemStart[sField]], array, matArray, ctx);
          CHKERRQ(ierr);
        }
#ifdef PETSC_USE_BOPT_g
        ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                     CHKERRQ(ierr);
#endif
      }
    }

    /* Setup global numbering, with reduction if necessary */
    ierr = GridCalcGeneralElementVecIndices(grid, elem, order, PETSC_NULL, PETSC_FALSE, vec);             CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
    ierr = PetscOptionsHasName(PETSC_NULL, "-trace_vec_assembly", &opt);                                  CHKERRQ(ierr);
    if (opt == PETSC_TRUE) {
      int var;
      for(var = 0; var < vec->reduceSize; var++)
        PetscPrintf(PETSC_COMM_SELF, "%2d %4.2g\n", vec->indices[var], PetscRealPart(array[var]));
    }
#endif
    /* Put values in global vector */
    ierr = ElementVecSetValues(vec, f, ADD_VALUES);                                                       CHKERRQ(ierr);
  }

  /* Evaluate self-interaction of new fields created by constraints */
  if (explicitConstraints == PETSC_TRUE) { 
    /* WARNING: This only accomodates 1 constrained field */
    /* Get constraint information */
    for(fieldIndex = 0; fieldIndex < numFields; fieldIndex++) {
      sField = fields[fieldIndex];
      if (grid->fields[sField].isConstrained == PETSC_TRUE) {
        newComp = grid->fields[sField].numComp + grid->fields[sField].constraintCompDiff;
        break;
      }
    }
    /* Calculate self-interaction */
    for(newField = 0; newField < numNewFields; newField++) {
      /* Initialize element matrix and vector */
      ierr = ElementMatZero(mat);                                                                         CHKERRQ(ierr);
      ierr = ElementVecZero(vec);                                                                         CHKERRQ(ierr);
      mat->reduceRowSize     = newComp;
      mat->reduceColSize     = newComp;
      elemAppVec->reduceSize = newComp;
      vec->reduceSize        = newComp;

      /* Calculate the indices and contribution to the element matrix from the new field */
      ierr = (*constCtx->ops->newelemmat)(constCtx, order, newField, mat);                                CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscOptionsHasName(PETSC_NULL, "-trace_mat_assembly_constrained", &opt);                    CHKERRQ(ierr);
      if (opt == PETSC_TRUE) {
        ierr = ElementMatView(mat, PETSC_VIEWER_STDOUT_(mat->comm));                                      CHKERRQ(ierr);
      }
#endif
      /* Global vector indices are the same as the matrix indices */
      for(row = 0; row < mat->reduceRowSize; row++) {
        vec->indices[row] = mat->rowIndices[row];
      }
      /* Local vector indices can be calculated directly from the field number */
      elemAppVec->indices[0] = grid->constraintOrder->firstVar[rank+1] - (numNewFields - newField)*newComp;
      for(row = 1; row < elemAppVec->reduceSize; row++) {
        elemAppVec->indices[row] = elemAppVec->indices[row-1]+1;
      }
      /* Retrieve element vector values from x */
      ierr = GridLocalToElementGeneral(grid, appVec, grid->bdReduceVecCur, reduceSystem, reduceElement, elemAppVec);CHKERRQ(ierr);
      /* Multiply element matrix and element vector */
      for(row = 0; row < mat->reduceRowSize; row++) {
        for(col = 0; col < mat->reduceColSize; col++) {
          vec->array[row] += mat->array[row*mat->reduceColSize+col]*elemAppVec->array[col];
        }
      }
      PetscLogFlops(2*mat->reduceRowSize*mat->reduceColSize);
      /* Put values in global vector */
      ierr = ElementVecSetValues(vec, f, ADD_VALUES);                                                     CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                        CHKERRQ(ierr);
#endif
    }
  }

  /* Reset element vector */
  elemGhostVec->reduceSize = locOrder->elemSize;

  ierr = VecDestroy(appVec);                                                                              CHKERRQ(ierr);
  ierr = ElementVecDestroy(elemAppVec);                                                                   CHKERRQ(ierr);
  if (grid->ALEActive == PETSC_TRUE) {
    ierr = ElementVecDestroy(ALEVec);                                                                     CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(f);                                                                             CHKERRQ(ierr);
  ierr = VecAssemblyEnd(f);                                                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GVecEvaluateSystemMatrixDiagonal_Triangular_2D"
int GVecEvaluateSystemMatrixDiagonal_Triangular_2D(Grid grid, GVec x, GVec d, void *ctx)
{
  Mesh                  mesh          = grid->mesh;
  int                   numNewFields  = grid->numNewFields;     /* The number of new fields added by constraints */
  int                   numMatOps     = grid->numMatOps;        /* The number of operators in the matrix */
  GridOp               *matOps        = grid->matOps;           /* The operators in the system matrix */
  VarOrdering           constOrder    = grid->constraintOrder;  /* The constrained variable ordering */
  PetscTruth            reduceSystem  = grid->reduceSystem;
  PetscTruth            reduceElement = grid->reduceElement;
  PetscTruth            expConst      = grid->explicitConstraints;
  PetscConstraintObject constCtx      = grid->constraintCtx;    /* The constraint object */
  int                   numFields     = grid->cm->numFields;    /* The number of fields in the calculation */
  int                  *fields        = grid->cm->fields;       /* The fields participating in the calculation */
  LocalVarOrdering      locOrder      = grid->locOrder;         /* The default local variable ordering */
  int                   elemSize      = locOrder->elemSize;     /* The number of shape functions in the element matrix */
  int                  *elemStart     = locOrder->elemStart;    /* The offset of each field in the element matrix */
  ElementMat            mat           = grid->mat;              /* The element matrix */
  PetscScalar          *array         = mat->array;             /* The values in the element matrix */
  Vec                   ghostVec      = grid->ghostVec;         /* The local solution vector */
  ElementVec            elemGhostVec  = grid->ghostElementVec;  /* The element vector from ghostVec */
  PetscScalar          *ghostArray    = elemGhostVec->array;    /* The values in elemGhostVec */
  MeshMover             mover;
  Grid                  ALEGrid;                                /* The grid describing the mesh velocity */
  VarOrdering           order;                                  /* The default variable ordering */
  ElementVec            MeshALEVec;                             /* ALE velocity vector with mesh discretization */
  ElementVec            ALEVec;                                 /* ALE velocity vector */
  PetscScalar          *ALEArray;                               /* The values in the ALE element vector */
  int                   newComp = 0;
  int                   numElements;
  int                   elem, f, sField, tField, op, newField;
#ifdef PETSC_USE_BOPT_g
  PetscTruth            opt;
#endif
  int                   ierr;

  PetscFunctionBegin;
  ierr = MeshGetMover(mesh, &mover);                                                                      CHKERRQ(ierr);
  ierr = MeshMoverGetVelocityGrid(mover, &ALEGrid);                                                       CHKERRQ(ierr);
  ierr = MeshGetInfo(mesh, PETSC_NULL, PETSC_NULL, PETSC_NULL, &numElements);                             CHKERRQ(ierr);
  if (expConst == PETSC_TRUE) {
    order = grid->constraintOrder;
  } else {
    order = grid->order;
  }
  /* Fill the local solution vectors */
  if (x != PETSC_NULL) {
    ierr = GridGlobalToLocal(grid, INSERT_VALUES, x);                                                     CHKERRQ(ierr);
  }

  /* Setup ALE variables -- No new variables should be ALE so ALEVec is not recalculated */
  if (grid->ALEActive == PETSC_TRUE) {
    /* Notice that the ALEArray is from this grid, not the mesh velocity grid */
    MeshALEVec = ALEGrid->vec;
    ierr       = ElementVecDuplicate(grid->vec, &ALEVec);                                                 CHKERRQ(ierr);
    ALEArray   = ALEVec->array;
  } else {
    MeshALEVec = PETSC_NULL;
    ALEArray   = PETSC_NULL;
  }

  /* Loop over elements */
  for(elem = 0; elem < numElements; elem++) {
    /* Initialize element matrix */
    ierr = ElementMatZero(mat);                                                                           CHKERRQ(ierr);
    mat->reduceRowSize       = locOrder->elemSize;
    mat->reduceColSize       = locOrder->elemSize;
    elemGhostVec->reduceSize = locOrder->elemSize;

    /* Setup local row indices for the ghost vector */
    ierr = GridCalcLocalElementVecIndices(grid, elem, elemGhostVec);                                      CHKERRQ(ierr);
    /* Setup local solution vector */
    ierr = GridLocalToElementGeneral(grid, ghostVec, grid->bdReduceVecCur, reduceSystem, reduceElement, elemGhostVec);CHKERRQ(ierr);
    /* Must transform to unconstrained variables for element integrals */
    ierr = GridProjectElementVec(grid, mesh, elem, order, PETSC_FALSE, elemGhostVec);                     CHKERRQ(ierr);

    /* Setup ALE variables */
    if (grid->ALEActive == PETSC_TRUE) {
      ierr = GridCalcLocalElementVecIndices(ALEGrid, elem, MeshALEVec);                                   CHKERRQ(ierr);
      ierr = GridLocalToElement(ALEGrid, MeshALEVec);                                                     CHKERRQ(ierr);
    }

    /* Calculate the contribution to the element matrix from each field */
    for(op = 0; op < numMatOps; op++) {
      sField = matOps[op].field;
      tField = grid->fields[sField].disc->operators[matOps[op].op]->test->field;
      if (grid->fields[sField].isActive) {
        if (matOps[op].isALE) {
          ierr = GridInterpolateElementVec(ALEGrid, 0, MeshALEVec, grid, sField, ALEVec);                 CHKERRQ(ierr);
          ierr = DiscretizationEvaluateALEOperatorGalerkin(grid->fields[sField].disc, mesh, elemSize, elemStart[tField], elemStart[sField],
                                                           matOps[op].op, matOps[op].alpha, elem, &ghostArray[elemStart[sField]],
                                                           ALEArray, array, ctx);
          CHKERRQ(ierr);
        } else {
          ierr = DiscretizationEvaluateOperatorGalerkin(grid->fields[sField].disc, mesh, elemSize, elemStart[tField], elemStart[sField],
                                                        matOps[op].op, matOps[op].alpha, elem, &ghostArray[elemStart[sField]],
                                                        array, ctx);
          CHKERRQ(ierr);
        }
#ifdef PETSC_USE_BOPT_g
        ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                      CHKERRQ(ierr);
#endif
      }
    }

    /* Setup global numbering, with reduction if necessary */
    ierr = GridCalcGeneralElementMatIndices(grid, elem, order, order, PETSC_FALSE, mat);                  CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
    ierr = PetscOptionsHasName(PETSC_NULL, "-trace_mat_assembly", &opt);                                  CHKERRQ(ierr);
    if (opt == PETSC_TRUE) {
      ierr = ElementMatView(mat, PETSC_VIEWER_STDOUT_(mat->comm));                                        CHKERRQ(ierr);
    }
#endif
    /* Put diagonal values in the global matrix */
    ierr = ElementMatSetDiagonalValues(mat, d, ADD_VALUES);                                               CHKERRQ(ierr);
  }

  /* Evaluate self-interaction of new fields created by constraints */
  if (expConst == PETSC_TRUE) { 
    /* WARNING: This only accomodates 1 constrained field */
    /* Get constraint information */
    for(f = 0; f < numFields; f++) {
      sField = fields[f];
      if (grid->fields[sField].isConstrained == PETSC_TRUE) {
        newComp = grid->fields[sField].numComp + grid->fields[sField].constraintCompDiff;
        break;
      }
    }
    /* Calculate self-interaction */
    for(newField = 0; newField < numNewFields; newField++) {
      /* Initialize element matrix */
      ierr = ElementMatZero(mat);                                                                         CHKERRQ(ierr);
      mat->reduceRowSize = newComp;
      mat->reduceColSize = newComp;

      /* Calculate the indices and contribution to the element matrix from the new field */
      ierr = (*constCtx->ops->newelemmat)(constCtx, constOrder, newField, mat);                           CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscOptionsHasName(PETSC_NULL, "-trace_mat_assembly_constrained", &opt);                    CHKERRQ(ierr);
      if (opt == PETSC_TRUE) {
        ierr = ElementMatView(mat, PETSC_VIEWER_STDOUT_(mat->comm));                                      CHKERRQ(ierr);
      }
#endif
      /* Put values in global matrix */
      ierr = ElementMatSetDiagonalValues(mat, d, ADD_VALUES);                                             CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
      ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                        CHKERRQ(ierr);
#endif
    }
  }

  /* Assemble matrix */
  ierr = VecAssemblyBegin(d);                                                                             CHKERRQ(ierr);
  ierr = VecAssemblyEnd(d);                                                                               CHKERRQ(ierr);

  /* Reset element matrix and vector */
  mat->reduceRowSize       = locOrder->elemSize;
  mat->reduceColSize       = locOrder->elemSize;
  elemGhostVec->reduceSize = locOrder->elemSize;

  /* Cleanup */
  if (grid->ALEActive == PETSC_TRUE) {
    ierr = ElementVecDestroy(ALEVec);                                                                     CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}