#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: fieldClassMap2d.c,v 1.6 2000/01/30 18:30:46 huangp Exp $";
#endif

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

/*---------------------------------------- Creation and Destruction Functions ---------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "FieldClassMapCreateTriangular2D"
/*@C
  FieldClassMapCreateTriangular2D - This function creates a class map for the given grid.

  Collective on Grid

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

  Output Parameter:
. map       - The map

  Level: beginner

.keywords: class, field class, class map, create
.seealso: FieldClassMapCreate(), FieldClassMapDestroy()
@*/
int FieldClassMapCreateTriangular2D(Grid grid, int numFields, int *fields, FieldClassMap *map)
{
  MPI_Comm      comm;
  ParameterDict dict;
  int           numTotalFields, f;
  int           ierr;

  PetscFunctionBegin;
  ierr = GridGetNumFields(grid, &numTotalFields);                                                         CHKERRQ(ierr);
  if (numFields > numTotalFields) {
    SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid number  %d of fields, grid only has %d", numFields, numTotalFields);
  }
  for(f = 0; f < numFields; f++) {
    if ((fields[f] < 0) || (fields[f] >= numTotalFields)) {
      SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", fields[f]);
    }
  }
  ierr = PetscObjectGetComm((PetscObject) grid, &comm);                                                   CHKERRQ(ierr);
  ierr = FieldClassMapCreate(comm, map);                                                                  CHKERRQ(ierr);
  ierr = PetscObjectGetComm((PetscObject) grid, &comm);                                                   CHKERRQ(ierr);
  ierr = ParameterDictCreate(comm, &dict);                                                                CHKERRQ(ierr);
  ierr = ParameterDictSetInteger(dict, "numFields", numFields);                                           CHKERRQ(ierr);
  ierr = ParameterDictSetObject(dict,  "fields",    fields);                                              CHKERRQ(ierr);
  ierr = ParameterDictSetObject(dict,  "grid",      grid);                                                CHKERRQ(ierr);
  ierr = PetscObjectSetParameterDict((PetscObject) *map, dict);                                           CHKERRQ(ierr);
  ierr = ParameterDictDestroy(dict);                                                                      CHKERRQ(ierr);
  ierr = FieldClassMapSetType(*map, CLASS_MAP_TRIANGULAR_2D);                                             CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "FieldClassMapDestroy_Triangular_2D"
int FieldClassMapDestroy_Triangular_2D(FieldClassMap map)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "FieldClassMapConstrain_Triangular_2D"
int FieldClassMapConstrain_Triangular_2D(FieldClassMap map, Grid grid, PetscTruth useBC, PetscTruth useConst,
                                         FieldClassMap *newMap)
{
  FieldClassMap cm;
  Mesh          mesh;
  Partition     part;
  int           numBd              = grid->numBd;
  int           numBC              = grid->numBC;
  GridBC       *gBC                = grid->bc;
  int           numPointBC         = grid->numPointBC;
  GridBC       *pBC                = grid->pointBC;
  int           numFields          = map->numFields;
  int          *fields             = map->fields;
  int           numNodes           = map->numNodes;
  int           numOverlapNodes    = map->numOverlapNodes;
  int           numClasses         = map->numClasses;
  int          *classes            = map->classes;
  int        ***localClassMap, ***classMap;
  PetscTruth    isConst;
  int           f, g, field, comp, node, bd, oldBd, bc, marker, nclass, newNClass, oldNClass;
  int           ierr;

  PetscFunctionBegin;
  ierr = FieldClassMapCreate(map->comm, &cm);                                                             CHKERRQ(ierr);
  ierr = PetscMemcpy(cm->ops, map->ops, sizeof(FieldClassMapOps));                                        CHKERRQ(ierr);
  ierr = PetscStrallocpy(map->serialize_name, &cm->serialize_name);                                       CHKERRQ(ierr);

  ierr = GridGetMesh(grid, &mesh);                                                                        CHKERRQ(ierr);
  cm->numFields       = numFields;
  cm->numNodes        = numNodes;
  cm->numOverlapNodes = numOverlapNodes;
  cm->numGhostNodes   = cm->numOverlapNodes - cm->numNodes;
  ierr = PetscMalloc(numFields * sizeof(int), &cm->fields);                                               CHKERRQ(ierr);
  ierr = PetscMemcpy(cm->fields, fields, numFields * sizeof(int));                                        CHKERRQ(ierr);
  /*
    classMap[numBd+numPointBC+1][numClasses][numFields]:

    The classMap[][][] structure is organized as follows. For every boundary, we determine the classes which are
    present on that boundary. For each class, we create a new class which does not contain the constrained fields.
    The classMap[][][] structure has for each boundary

      flag_0, flag_1, \ldots, flag_numFields-1
      flag_0, flag_1, \ldots, flag_numFields-1
      ...
      flag_0, flag_1, \ldots, flag_numFields-1

    where a positive flag indicates that the field is constrained in that class, and a zero indicates no change.
    Point boundary conditions are treated like additional boundaries, and since we only allow one constraint,
    the last position is used for the constraint.
  */
  cm->mapSize       = numBd + numPointBC + 1;
  cm->numOldClasses = numClasses;
  ierr = PetscMalloc(cm->mapSize * sizeof(int **), &localClassMap);                                       CHKERRQ(ierr);
  ierr = PetscMalloc(cm->mapSize * sizeof(int **), &classMap);                                            CHKERRQ(ierr);
  for(bd = 0; bd < cm->mapSize; bd++) {
    ierr = PetscMalloc(numClasses * sizeof(int *), &localClassMap[bd]);                                   CHKERRQ(ierr);
    ierr = PetscMalloc(numClasses * sizeof(int *), &classMap[bd]);                                        CHKERRQ(ierr);
    for(nclass = 0; nclass < numClasses; nclass++) {
      ierr = PetscMalloc(numFields * sizeof(int), &localClassMap[bd][nclass]);                            CHKERRQ(ierr);
      ierr = PetscMalloc(numFields * sizeof(int), &classMap[bd][nclass]);                                 CHKERRQ(ierr);
      ierr = PetscMemzero(localClassMap[bd][nclass], numFields * sizeof(int));                            CHKERRQ(ierr);
      ierr = PetscMemzero(classMap[bd][nclass],      numFields * sizeof(int));                            CHKERRQ(ierr);
    }
  }
  /* Mark constrained classes for reassignment */
  for(node = 0; node < numNodes; node++) {
    ierr = MeshGetNodeBoundary(mesh, node, &marker);                                                      CHKERRQ(ierr);
    if (marker == 0)
      continue;
    nclass = classes[node];
    if ((useConst == PETSC_TRUE) && (marker < 0)) {
      /* All constrained nodes have negative markers */
      for(f = 0; f < numFields; f++) {
        field = fields[f];
        if (grid->fields[field].isConstrained == PETSC_TRUE) {
          localClassMap[cm->mapSize-1][nclass][f] = 1;
        }
      }
    } else if ((useBC == PETSC_TRUE) && (marker > 0)) {
      for(bc = 0; bc < numBC; bc++) {
        if (gBC[bc].reduce == PETSC_FALSE) continue;
        if (gBC[bc].boundary == marker) {
          ierr = MeshGetBoundaryIndex(mesh, marker, &bd);                                                 CHKERRQ(ierr);
          for(f = 0; f < numFields; f++) {
            if (fields[f] == gBC[bc].field) {
              localClassMap[bd][nclass][f] = 1;
              break;
            }
          }
        }
      }
      for(bc = 0; bc < numPointBC; bc++) {
        if (pBC[bc].reduce == PETSC_FALSE) continue;
        if (pBC[bc].node == node) {
          for(f = 0; f < numFields; f++) {
            if (fields[f] == pBC[bc].field) {
              localClassMap[numBd+bc][nclass][f] = 1;
              break;
            }
          }
        }
      }
    }
  }
  /* Synchronize localClassMap */
  for(bd = 0; bd < cm->mapSize; bd++) {
    for(nclass = 0; nclass < numClasses; nclass++) {
      ierr = MPI_Allreduce(localClassMap[bd][nclass], classMap[bd][nclass], numFields, MPI_INT, MPI_LOR, cm->comm);
      CHKERRQ(ierr);
    }
  }
  /* Assign new classes */
  ierr = PetscMalloc(cm->mapSize * sizeof(int *), &cm->classMap);                                         CHKERRQ(ierr);
  for(bd = 0; bd < cm->mapSize; bd++) {
    ierr = PetscMalloc(cm->numOldClasses * sizeof(int), &cm->classMap[bd]);                               CHKERRQ(ierr);
    for(nclass = 0; nclass < numClasses; nclass++) {
      cm->classMap[bd][nclass] = nclass;
    }
  }
  for(bd = 0, cm->numClasses = numClasses; bd < cm->mapSize; bd++) {
    for(nclass = 0; nclass < numClasses; nclass++) {
      for(f = 0, isConst = PETSC_FALSE; f < numFields; f++) {
        if (classMap[bd][nclass][f] != 0) {
          isConst = PETSC_TRUE;
          field   = fields[f];
          /* Map old classes to new with the condition which created it for a given field */
          classMap[bd][nclass][f] = cm->numClasses;
          /* Check that field should be constrained */
          if ((bd == cm->mapSize-1) && (grid->fields[field].isConstrained == PETSC_FALSE)) {
            SETERRQ1(PETSC_ERR_PLIB, "Invalid constrained field %d", field);
          }
          /* Must include previous constraints on the node */
          for(g = 0; g < numFields; g++) {
            for(oldBd = 0; oldBd < bd; oldBd++) {
              for(oldNClass = 0; oldNClass < numClasses; oldNClass++) {
                if (classMap[oldBd][oldNClass][g] != 0) {
                  classMap[bd][nclass][g] = cm->numClasses;
                  break;
                }
              }
            }
          }
          break;
        }
      }
      /* Create a new class */
      if (isConst == PETSC_TRUE) {
        cm->classMap[bd][nclass] = cm->numClasses++;
      }
    }
  }
  /* Set new classes on nodes */
  ierr = PetscMalloc(numOverlapNodes * sizeof(int), &cm->classes);                                        CHKERRQ(ierr);
  ierr = PetscMemcpy(cm->classes, classes, numOverlapNodes * sizeof(int));                                CHKERRQ(ierr);
  for(node = 0; node < numOverlapNodes; node++) {
    ierr = MeshGetNodeBoundary(mesh, node, &marker);                                                      CHKERRQ(ierr);
    if (marker == 0)
      continue;
    nclass = classes[node];
    if (marker < 0) {
      cm->classes[node] = cm->classMap[cm->mapSize-1][nclass];
    } else {
      for(bc = 0; bc < numBC; bc++) {
        if (gBC[bc].reduce == PETSC_FALSE) continue;
        if (gBC[bc].boundary == marker) {
          ierr = MeshGetBoundaryIndex(mesh, marker, &bd);                                                 CHKERRQ(ierr);
          if (cm->classMap[bd][nclass] != nclass) {
            cm->classes[node] = cm->classMap[bd][nclass];
          }
        }
      }
      /* Notice that point BC must come last since these nodes may lie on a boundary which was also constrained */
      for(bc = 0; bc < numPointBC; bc++) {
        if (pBC[bc].reduce == PETSC_FALSE) continue;
        if (pBC[bc].node == node) {
          if (cm->classMap[numBd+bc][nclass] != nclass) {
            cm->classes[node] = cm->classMap[numBd+bc][nclass];
          }
        }
      }
    }
  }
  /* Calculate new class structure */
  cm->isReduced          = map->isReduced;
  cm->isConstrained      = map->isConstrained;
  ierr = PetscMalloc(numFields      * sizeof(int *), &cm->fieldClasses);                                  CHKERRQ(ierr);
  ierr = PetscMalloc(numFields      * sizeof(int *), &cm->reduceFieldClasses);                            CHKERRQ(ierr);
  ierr = PetscMalloc(cm->numClasses * sizeof(int),   &cm->classSizes);                                    CHKERRQ(ierr);
  ierr = PetscMalloc(cm->numClasses * sizeof(int),   &cm->isClassConstrained);                            CHKERRQ(ierr);
  ierr = PetscMalloc(cm->numClasses * sizeof(int),   &cm->classSizeDiffs);                                CHKERRQ(ierr);
  ierr = PetscMemzero(cm->isClassConstrained, cm->numClasses * sizeof(int));                              CHKERRQ(ierr);
  ierr = PetscMemzero(cm->classSizeDiffs,     cm->numClasses * sizeof(int));                              CHKERRQ(ierr);
  /* Copy old structure */
  for(nclass = 0; nclass < numClasses; nclass++) {
    cm->classSizes[nclass]         = map->classSizes[nclass];
    cm->isClassConstrained[nclass] = map->isClassConstrained[nclass];
    cm->classSizeDiffs[nclass]     = map->classSizeDiffs[nclass];
  }
  for(newNClass = numClasses; newNClass < cm->numClasses; newNClass++) {
    /* Invert the class map */
    for(bd = 0; bd < cm->mapSize; bd++) {
      for(nclass = 0; nclass < numClasses; nclass++)
        if (cm->classMap[bd][nclass] == newNClass)
          break;
      if (nclass < numClasses)
        break;
    }
    if ((nclass >= numClasses) || (bd >= cm->mapSize)) {
      SETERRQ1(PETSC_ERR_PLIB, "Unable to locate class %d in class map", newNClass);
    }

    cm->classSizes[newNClass]     = map->classSizes[nclass];
    cm->classSizeDiffs[newNClass] = map->classSizeDiffs[nclass];
  }
  /* Apply constraints */
  for(f = 0; f < numFields; f++) {
    ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);                               CHKERRQ(ierr);
    ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->reduceFieldClasses[f]);                         CHKERRQ(ierr);
    ierr = PetscMemzero(cm->reduceFieldClasses[f], cm->numClasses * sizeof(int));                         CHKERRQ(ierr);
    /* Copy old structure */
    for(nclass = 0; nclass < numClasses; nclass++) {
      cm->fieldClasses[f][nclass]         = map->fieldClasses[f][nclass];
      if (map->reduceFieldClasses) {
        cm->reduceFieldClasses[f][nclass] = map->reduceFieldClasses[f][nclass];
      }
    }

    /* Remove constrained fields from new classes */
    for(newNClass = numClasses; newNClass < cm->numClasses; newNClass++) {
      /* Invert the class map */
      for(bd = 0; bd < cm->mapSize; bd++) {
        for(nclass = 0; nclass < numClasses; nclass++)
          if (cm->classMap[bd][nclass] == newNClass)
            break;
        if (nclass < numClasses)
          break;
      }
      if ((nclass >= numClasses) || (bd >= cm->mapSize)) {
        SETERRQ1(PETSC_ERR_PLIB, "Unable to locate class %d in class map", newNClass);
      }
      /* Check every BC for constraints on this field/class */
      cm->fieldClasses[f][newNClass]             = map->fieldClasses[f][nclass];
      for(bd = 0; bd < cm->mapSize; bd++) {
        if (classMap[bd][nclass][f] == newNClass) {
          /* New classes can be redefined by another BC to be a higher number, but do not
             operate on the existing new class with that BC (hence the second check) */
          if (cm->classMap[bd][nclass] != newNClass) {
            SETERRQ3(PETSC_ERR_PLIB, "Invalid bc %d constraining class %d to %d", bd, nclass, newNClass);
          }
          /* If the field is constrained, it has no values on these nodes */
          cm->fieldClasses[f][newNClass]         = 0;
          if (bd < cm->mapSize-1) {
            cm->isReduced                        = PETSC_TRUE;
            cm->reduceFieldClasses[f][newNClass] = 1;
          }
          /* Subtract the degrees of freedom of the constrained field for each BC */
          ierr = GridGetFieldComponents(grid, fields[f], &comp);                                          CHKERRQ(ierr);
          cm->classSizes[newNClass]             -= comp;
          if (cm->classSizes[newNClass] < 0) {
            SETERRQ2(PETSC_ERR_PLIB, "Invalid size for class %d (derived from class %d)", newNClass, nclass);
          }
          /* Handle user constraints separately from BC */
          if (bd == cm->mapSize-1) {
            cm->isConstrained                    = PETSC_TRUE;
            /* This class is constrained */
            cm->isClassConstrained[newNClass]    = 1;
            /* There are this many components in the new field */
            cm->classSizeDiffs[newNClass]       += grid->fields[fields[f]].constraintCompDiff + comp;
          }
        }
      }
    }
  }

  /* Cleanup */
  if (cm->isReduced == PETSC_FALSE) {
    for(f = 0; f < map->numFields; f++) {
      ierr = PetscFree(cm->reduceFieldClasses[f]);                                                        CHKERRQ(ierr);
    }
    ierr = PetscFree(cm->reduceFieldClasses);                                                             CHKERRQ(ierr);
  }
  for(bd = 0; bd < cm->mapSize; bd++) {
    for(nclass = 0; nclass < numClasses; nclass++) {
      ierr = PetscFree(localClassMap[bd][nclass]);                                                        CHKERRQ(ierr);
      ierr = PetscFree(classMap[bd][nclass]);                                                             CHKERRQ(ierr);
    }
    ierr = PetscFree(localClassMap[bd]);                                                                  CHKERRQ(ierr);
    ierr = PetscFree(classMap[bd]);                                                                       CHKERRQ(ierr);
  }
  ierr = PetscFree(localClassMap);                                                                        CHKERRQ(ierr);
  ierr = PetscFree(classMap);                                                                             CHKERRQ(ierr);

  /* Communicate ghost classes */
  ierr = MeshGetPartition(mesh, &part);                                                                   CHKERRQ(ierr);
  ierr = PartitionGhostNodeExchange(part, INSERT_VALUES, SCATTER_FORWARD, cm->classes, &cm->classes[cm->numNodes]);
  CHKERRQ(ierr);

  PetscLogObjectMemory(cm, (cm->mapSize + numFields*2) * sizeof(int *) + (numFields + cm->numOldClasses*cm->mapSize +
                       cm->numClasses*(numFields*2 + 3) + numOverlapNodes) * sizeof(int));
  *newMap = cm;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "FieldClassMapReduce_Triangular_2D"
int FieldClassMapReduce_Triangular_2D(FieldClassMap map, Grid grid, FieldClassMap *newMap) {
  FieldClassMap cm;
  Mesh          mesh;
  Partition     part;
  int           numBd              = grid->numBd;
  int           numBC              = grid->numBC;
  GridBC       *gBC                = grid->bc;
  int           numPointBC         = grid->numPointBC;
  GridBC       *pBC                = grid->pointBC;
  int           numFields          = map->numFields;
  int          *fields             = map->fields;
  int           numNodes           = map->numNodes;
  int           numOverlapNodes    = map->numOverlapNodes;
  int           numClasses         = map->numClasses;
  int          *classes            = map->classes;
  int        ***localClassMap, ***classMap;
  PetscTruth    isConst;
  int           f, g, field, comp, node, bd, bc, marker, nclass, newNClass;
  int           ierr;

  PetscFunctionBegin;
  ierr = FieldClassMapCreate(map->comm, &cm);                                                             CHKERRQ(ierr);
  ierr = PetscMemcpy(cm->ops, map->ops, sizeof(FieldClassMapOps));                                        CHKERRQ(ierr);
  ierr = PetscStrallocpy(map->serialize_name, &cm->serialize_name);                                       CHKERRQ(ierr);

  ierr = GridGetMesh(grid, &mesh);                                                                        CHKERRQ(ierr);
  cm->numNodes        = numNodes;
  cm->numOverlapNodes = numOverlapNodes;
  cm->numGhostNodes   = cm->numOverlapNodes - cm->numNodes;
  ierr = PetscMalloc(numFields * sizeof(int), &cm->fields);                                               CHKERRQ(ierr);
  /*
    classMap[numBd+numPointBC][numClasses][numFields]:

    The classMap[][][] structure is organized as follows. For every boundary, we determine the classes which are
    present on that boundary. For each class, we create a new class which does not contain the constrained fields.
    The classMap[][][] structure has for each boundary

      flag_0, flag_1, \ldots, flag_numFields-1
      flag_0, flag_1, \ldots, flag_numFields-1
      ...
      flag_0, flag_1, \ldots, flag_numFields-1

    where a positive flag indicates that the field is constrained in that class, and a zero indicates no change.
    Point boundary conditions are treated like additional boundaries.
  */
  cm->mapSize       = numBd + numPointBC;
  cm->numOldClasses = numClasses;
  ierr = PetscMalloc(cm->mapSize * sizeof(int **), &localClassMap);                                       CHKERRQ(ierr);
  ierr = PetscMalloc(cm->mapSize * sizeof(int **), &classMap);                                            CHKERRQ(ierr);
  for(bd = 0; bd < cm->mapSize; bd++) {
    ierr = PetscMalloc(numClasses * sizeof(int *), &localClassMap[bd]);                                   CHKERRQ(ierr);
    ierr = PetscMalloc(numClasses * sizeof(int *), &classMap[bd]);                                        CHKERRQ(ierr);
    for(nclass = 0; nclass < numClasses; nclass++) {
      ierr = PetscMalloc(numFields * sizeof(int), &localClassMap[bd][nclass]);                            CHKERRQ(ierr);
      ierr = PetscMalloc(numFields * sizeof(int), &classMap[bd][nclass]);                                 CHKERRQ(ierr);
      ierr = PetscMemzero(localClassMap[bd][nclass], numFields * sizeof(int));                            CHKERRQ(ierr);
      ierr = PetscMemzero(classMap[bd][nclass],      numFields * sizeof(int));                            CHKERRQ(ierr);
    }
  }
  /* Mark constrained classes for reassignment */
  for(node = 0; node < numNodes; node++) {
    ierr = MeshGetNodeBoundary(mesh, node, &marker);                                                      CHKERRQ(ierr);
    if (marker == 0)
      continue;
    nclass = classes[node];
    if (marker > 0) {
      for(bc = 0; bc < numBC; bc++) {
        if (gBC[bc].reduce == PETSC_FALSE)
          continue;
        if (gBC[bc].boundary == marker) {
          ierr = MeshGetBoundaryIndex(mesh, marker, &bd);                                                 CHKERRQ(ierr);
          for(f = 0; f < numFields; f++) {
            if (fields[f] == gBC[bc].field) {
              localClassMap[bd][nclass][f] = 1;
              break;
            }
          }
        }
      }
      for(bc = 0; bc < numPointBC; bc++) {
        if (pBC[bc].reduce == PETSC_FALSE) continue;
        if (pBC[bc].node == node) {
          for(f = 0; f < numFields; f++) {
            if (fields[f] == pBC[bc].field) {
              localClassMap[numBd+bc][nclass][f] = 1;
              break;
            }
          }
        }
      }
    }
  }
  /* Synchronize localClassMap */
  for(bd = 0; bd < cm->mapSize; bd++) {
    for(nclass = 0; nclass < numClasses; nclass++) {
      ierr = MPI_Allreduce(localClassMap[bd][nclass], classMap[bd][nclass], numFields, MPI_INT, MPI_LOR, cm->comm);
      CHKERRQ(ierr);
    }
  }

  /* Assign new classes and build field list */
  ierr = PetscMalloc(cm->mapSize * sizeof(int *), &cm->classMap);                                         CHKERRQ(ierr);
  for(bd = 0; bd < cm->mapSize; bd++) {
    ierr = PetscMalloc(cm->numOldClasses * sizeof(int), &cm->classMap[bd]);                               CHKERRQ(ierr);
    for(nclass = 0; nclass < numClasses; nclass++) {
      cm->classMap[bd][nclass] = nclass;
    }
  }
  cm->numFields = 0;
  for(bd = 0, cm->numClasses = numClasses; bd < cm->mapSize; bd++) {
    for(nclass = 0; nclass < numClasses; nclass++) {
      for(f = 0, isConst = PETSC_FALSE; f < numFields; f++) {
        if (classMap[bd][nclass][f] != 0) {
          isConst = PETSC_TRUE;
          field   = fields[f];
          /* Add to list of reduced fields */
          for(g = 0; g < cm->numFields; g++)
            if (field == cm->fields[g])
              break;
          if (g == cm->numFields)
            cm->fields[cm->numFields++] = field;
        }
      }
      /* Create a new class */
      if (isConst == PETSC_TRUE) {
        cm->classMap[bd][nclass] = cm->numClasses++;
      }
    }
  }

  /* Set new classes on nodes */
  ierr = PetscMalloc(numOverlapNodes * sizeof(int), &cm->classes);                                        CHKERRQ(ierr);
  ierr = PetscMemcpy(cm->classes, classes, numOverlapNodes * sizeof(int));                                CHKERRQ(ierr);
  for(node = 0; node < numOverlapNodes; node++) {
    ierr = MeshGetNodeBoundary(mesh, node, &marker);                                                      CHKERRQ(ierr);
    if (marker == 0)
      continue;
    nclass = classes[node];
    for(bc = 0; bc < numBC; bc++) {
      if (gBC[bc].reduce == PETSC_FALSE) continue;
      if (gBC[bc].boundary == marker) {
        ierr = MeshGetBoundaryIndex(mesh, marker, &bd);                                                   CHKERRQ(ierr);
        cm->classes[node] = cm->classMap[bd][nclass];
      }
    }
    /* Notice that point BC must come last since these nodes may lie on a boundary which was also constrained */
    for(bc = 0; bc < numPointBC; bc++) {
      if (pBC[bc].reduce == PETSC_FALSE) continue;
      if (pBC[bc].node == node) {
        cm->classes[node] = cm->classMap[numBd+bc][nclass];
      }
    }
  }
  /* Calculate new class structure */
  cm->isReduced          = PETSC_TRUE;
  cm->isConstrained      = map->isConstrained;
  ierr = PetscMalloc(cm->numFields  * sizeof(int *), &cm->fieldClasses);                                  CHKERRQ(ierr);
  ierr = PetscMalloc(cm->numFields  * sizeof(int *), &cm->reduceFieldClasses);                            CHKERRQ(ierr);
  ierr = PetscMalloc(cm->numClasses * sizeof(int),   &cm->isClassConstrained);                            CHKERRQ(ierr);
  ierr = PetscMalloc(cm->numClasses * sizeof(int),   &cm->classSizeDiffs);                                CHKERRQ(ierr);
  ierr = PetscMemzero(cm->isClassConstrained, cm->numClasses * sizeof(int));                              CHKERRQ(ierr);
  ierr = PetscMemzero(cm->classSizeDiffs,     cm->numClasses * sizeof(int));                              CHKERRQ(ierr);
  /* Copy old structure */
  for(nclass = 0; nclass < numClasses; nclass++) {
    cm->isClassConstrained[nclass] = map->isClassConstrained[nclass];
    cm->classSizeDiffs[nclass]     = map->classSizeDiffs[nclass];
  }
  /* Apply constraints */
  for(f = 0; f < cm->numFields; f++) {
    ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);                               CHKERRQ(ierr);
    ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->reduceFieldClasses[f]);                         CHKERRQ(ierr);
    ierr = PetscMemzero(cm->fieldClasses[f], cm->numClasses * sizeof(int));                               CHKERRQ(ierr);
    /* Copy old structure */
    for(nclass = 0; nclass < numClasses; nclass++) {
      cm->reduceFieldClasses[f][nclass] = map->fieldClasses[f][nclass];
      if (map->reduceFieldClasses) {
        cm->fieldClasses[f][nclass]     = map->reduceFieldClasses[f][nclass];
      }
    }

    /* Remove constrained fields from new classes */
    for(newNClass = numClasses; newNClass < cm->numClasses; newNClass++) {
      /* Invert the class map */
      for(bd = 0; bd < cm->mapSize; bd++) {
        for(nclass = 0; nclass < numClasses; nclass++)
          if (cm->classMap[bd][nclass] == newNClass)
            break;
        if (nclass < numClasses)
          break;
      }
      if ((nclass >= numClasses) || (bd >= cm->mapSize)) {
        SETERRQ1(PETSC_ERR_PLIB, "Unable to locate class %d in class map", newNClass);
      }
      /* Check every BC for constraints on this field/class */
      cm->reduceFieldClasses[f][newNClass]     = map->fieldClasses[f][nclass];
      for(bd = 0; bd < cm->mapSize; bd++) {
        /* New classes can be redefined by another BC to be a higher number, but do not
           operate on the existing new class with that BC (hence the second check)
        */
        if ((classMap[bd][nclass][f] == 1) && (cm->classMap[bd][nclass] <= newNClass)) {
          /* If the field is constrained, it has no values on these nodes */
          cm->fieldClasses[f][newNClass]       = 1;
          cm->reduceFieldClasses[f][newNClass] = 0;
        }
      }
    }
  }

  /* Setup class sizes */
  ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizes);                                      CHKERRQ(ierr);
  ierr = PetscMemzero(cm->classSizes, cm->numClasses * sizeof(int));                                      CHKERRQ(ierr);
  for(nclass = 0; nclass < cm->numClasses; nclass++) {
    for(f = 0; f < cm->numFields; f++) {
      field = cm->fields[f];
      if (cm->fieldClasses[f][nclass]) {
        ierr = GridGetFieldComponents(grid, field, &comp);                                                CHKERRQ(ierr);
        cm->classSizes[nclass] += comp;
      }
    }
  }

  /* Cleanup */
  for(bd = 0; bd < cm->mapSize; bd++) {
    for(nclass = 0; nclass < numClasses; nclass++) {
      ierr = PetscFree(localClassMap[bd][nclass]);                                                        CHKERRQ(ierr);
      ierr = PetscFree(classMap[bd][nclass]);                                                             CHKERRQ(ierr);
    }
    ierr = PetscFree(localClassMap[bd]);                                                                  CHKERRQ(ierr);
    ierr = PetscFree(classMap[bd]);                                                                       CHKERRQ(ierr);
  }
  ierr = PetscFree(localClassMap);                                                                        CHKERRQ(ierr);
  ierr = PetscFree(classMap);                                                                             CHKERRQ(ierr);

  /* Communicate ghost classes */
  ierr = MeshGetPartition(mesh, &part);                                                                   CHKERRQ(ierr);
  ierr = PartitionGhostNodeExchange(part, INSERT_VALUES, SCATTER_FORWARD, cm->classes, &cm->classes[cm->numNodes]);
  CHKERRQ(ierr);

  PetscLogObjectMemory(cm, (cm->mapSize + cm->numFields*2) * sizeof(int *) + (numFields + cm->numOldClasses*cm->mapSize +
                       cm->numClasses*(cm->numFields*2 + 3) + numOverlapNodes) * sizeof(int));
  *newMap = cm;
  PetscFunctionReturn(0);
}

static FieldClassMapOps cmOps = {0,
                                 FieldClassMapConstrain_Triangular_2D,
                                 FieldClassMapReduce_Triangular_2D,
                                 FieldClassMapDestroy_Triangular_2D,
                                 0};

EXTERN_C_BEGIN
#undef  __FUNCT__
#define __FUNCT__ "FieldClassMapSerialize_Triangular_2D"
int FieldClassMapSerialize_Triangular_2D(MPI_Comm comm, FieldClassMap *map, PetscViewer viewer, PetscTruth store)
{
  FieldClassMap cm;
  int           fd;
  int           zero = 0;
  int           one  = 0;
  int           hasClassMap;
  int           bd, f;
  int           ierr;

  PetscFunctionBegin;
  ierr = PetscViewerBinaryGetDescriptor(viewer, &fd);                                                     CHKERRQ(ierr);
  if (store) {
    cm   = *map;
    ierr = PetscBinaryWrite(fd, &cm->numFields,          1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  cm->fields,             cm->numFields,        PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &cm->numNodes,           1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &cm->numGhostNodes,      1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &cm->numOverlapNodes,    1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &cm->numClasses,         1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    for(f = 0; f < cm->numFields; f++) {
      ierr = PetscBinaryWrite(fd,  cm->fieldClasses[f],  cm->numClasses,       PETSC_INT, 0);             CHKERRQ(ierr);
    }
    ierr = PetscBinaryWrite(fd,  cm->classes,            cm->numOverlapNodes,  PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  cm->classSizes,         cm->numClasses,       PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &cm->isReduced,          1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    if (cm->isReduced == PETSC_TRUE) {
      for(f = 0; f < cm->numFields; f++) {
        ierr = PetscBinaryWrite(fd,  cm->reduceFieldClasses[f], cm->numClasses, PETSC_INT, 0);            CHKERRQ(ierr);
      }
    }
    ierr = PetscBinaryWrite(fd, &cm->isConstrained,      1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  cm->isClassConstrained, cm->numClasses,       PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  cm->classSizeDiffs,     cm->numClasses,       PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &cm->mapSize,            1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &cm->numOldClasses,      1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    if (cm->classMap != PETSC_NULL) {
      ierr = PetscBinaryWrite(fd, &one,                  1,                    PETSC_INT, 0);             CHKERRQ(ierr);
      for(bd = 0; bd < cm->mapSize; bd++) {
        ierr = PetscBinaryWrite(fd,  cm->classMap[bd],   cm->numOldClasses,    PETSC_INT, 0);             CHKERRQ(ierr);
      }
    } else {
      ierr = PetscBinaryWrite(fd, &zero,                 1,                    PETSC_INT, 0);             CHKERRQ(ierr);
    }
  } else {
    /* Create the mesh context */
    ierr = FieldClassMapCreate(comm, &cm);                                                                CHKERRQ(ierr);
    ierr = PetscMemcpy(cm->ops, &cmOps, sizeof(FieldClassMapOps));                                        CHKERRQ(ierr);
    ierr = PetscStrallocpy(CLASS_MAP_TRIANGULAR_2D, &cm->type_name);                                      CHKERRQ(ierr);
    cm->data = PETSC_NULL;

    ierr = PetscBinaryRead(fd, &cm->numFields,          1,                   PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscMalloc(cm->numFields * sizeof(int), &cm->fields);                                         CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd,  cm->fields,             cm->numFields,       PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &cm->numNodes,           1,                   PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &cm->numGhostNodes,      1,                   PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &cm->numOverlapNodes,    1,                   PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &cm->numClasses,         1,                   PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscMalloc(cm->numFields * sizeof(int *), &cm->fieldClasses);                                 CHKERRQ(ierr);
    for(f = 0; f < cm->numFields; f++) {
      ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);                             CHKERRQ(ierr);
      ierr = PetscBinaryRead(fd,  cm->fieldClasses[f],  cm->numClasses,      PETSC_INT);                  CHKERRQ(ierr);
    }
    ierr = PetscMalloc(cm->numOverlapNodes * sizeof(int), &cm->classes);                                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd,  cm->classes,            cm->numOverlapNodes, PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscMalloc(cm->numClasses      * sizeof(int), &cm->classSizes);CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd,  cm->classSizes,         cm->numClasses,      PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &cm->isReduced,          1,                   PETSC_INT);                  CHKERRQ(ierr);
    if (cm->isReduced == PETSC_TRUE) {
      ierr = PetscMalloc(cm->numFields * sizeof(int *), &cm->reduceFieldClasses);                         CHKERRQ(ierr);
      for(f = 0; f < cm->numFields; f++) {
        ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->reduceFieldClasses[f]);                     CHKERRQ(ierr);
        ierr = PetscBinaryRead(fd, cm->reduceFieldClasses[f], cm->numClasses, PETSC_INT);                 CHKERRQ(ierr);
      }
    }
    ierr = PetscBinaryRead(fd, &cm->isConstrained,      1,                   PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->isClassConstrained);                            CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd,  cm->isClassConstrained, cm->numClasses,      PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizeDiffs);                                CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd,  cm->classSizeDiffs,     cm->numClasses,      PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &cm->mapSize,            1,                   PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &cm->numOldClasses,      1,                   PETSC_INT);                  CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &hasClassMap,            1,                   PETSC_INT);                  CHKERRQ(ierr);
    if (hasClassMap) {
      ierr = PetscMalloc(cm->mapSize * sizeof(int *), &cm->classMap);                                     CHKERRQ(ierr);
      for(bd = 0; bd < cm->mapSize; bd++) {
        ierr = PetscMalloc(cm->numOldClasses * sizeof(int), &cm->classMap[bd]);                           CHKERRQ(ierr);
        ierr = PetscBinaryRead(fd,  cm->classMap[bd],   cm->numOldClasses,   PETSC_INT);                  CHKERRQ(ierr);
      }
    }

    PetscLogObjectMemory(cm, (cm->numFields + cm->mapSize) * sizeof(int *) + (cm->numFields*(cm->numClasses + 1) +
                         cm->numOverlapNodes + cm->numOldClasses*cm->mapSize + cm->numClasses*3) * sizeof(int));
    *map = cm;
  }
  PetscFunctionReturn(0);
}
EXTERN_C_END

EXTERN_C_BEGIN
#undef  __FUNCT__
#define __FUNCT__ "FieldClassMapCreate_Triangular_2D"
int FieldClassMapCreate_Triangular_2D(FieldClassMap cm)
{
  ParameterDict  dict;
  Grid           grid;
  Mesh           mesh;
  Partition      part;
  Discretization disc, firstDisc;
  int           *fields;
  int            numOverlapElements, numCorners;
  int            f, field, comp, elem, corner, node, nclass;
  PetscTruth     match, isconst, islin, isquad;
  int            ierr;

  PetscFunctionBegin;
  ierr = PetscMemcpy(cm->ops, &cmOps, sizeof(FieldClassMapOps));                                          CHKERRQ(ierr);
  ierr = PetscStrallocpy(CLASS_MAP_SER_TRIANGULAR_2D_BINARY, &cm->serialize_name);                        CHKERRQ(ierr);
  /* Get arguments */
  ierr = PetscObjectGetParameterDict((PetscObject) cm, &dict);                                            CHKERRQ(ierr);
  ierr = ParameterDictGetInteger(dict, "numFields", &cm->numFields);                                      CHKERRQ(ierr);
  ierr = ParameterDictGetObject(dict,  "fields",    (void **) &fields);                                   CHKERRQ(ierr);
  ierr = ParameterDictGetObject(dict,  "grid",      (void **) &grid);                                     CHKERRQ(ierr);
  /* Get fields */
  ierr = PetscMalloc(cm->numFields * sizeof(int), &cm->fields);                                           CHKERRQ(ierr);
  ierr = PetscMemcpy(cm->fields, fields, cm->numFields * sizeof(int));                                    CHKERRQ(ierr);

  /* Get number of classes -- this is very simplistic, needs to be changed for other discretizations */
  ierr = GridGetFieldDisc(grid, 0, &firstDisc);                                                           CHKERRQ(ierr);
  cm->numClasses = 1;
  for(f = 0; f < cm->numFields; f++) {
    ierr = GridGetFieldDisc(grid, cm->fields[f], &disc);                                                  CHKERRQ(ierr);
    ierr = PetscStrcmp(disc->type_name, firstDisc->type_name, &match);                                    CHKERRQ(ierr);
    if (match == PETSC_FALSE) {
      cm->numClasses++;
      break;
    }
  }
  ierr = GridGetMesh(grid, &mesh);                                                                        CHKERRQ(ierr);
  ierr = MeshGetPartition(mesh, &part);                                                                   CHKERRQ(ierr);
  ierr = PartitionGetNumNodes(part, &cm->numNodes);                                                       CHKERRQ(ierr);
  ierr = PartitionGetNumOverlapNodes(part, &cm->numOverlapNodes);                                         CHKERRQ(ierr);
  cm->numGhostNodes = cm->numOverlapNodes - cm->numNodes;
  ierr = PetscMalloc(cm->numOverlapNodes * sizeof(int), &cm->classes);                                    CHKERRQ(ierr);

  /* Setup field class membership -- this is very simplistic, needs to be changed for other discretizations */
  ierr = PetscMalloc(cm->numFields * sizeof(int *), &cm->fieldClasses);                                   CHKERRQ(ierr);
  for(f = 0; f < cm->numFields; f++) {
    field = cm->fields[f];
    ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->fieldClasses[f]);                               CHKERRQ(ierr);
    ierr = GridGetFieldDisc(grid, field, &disc);                                                          CHKERRQ(ierr);
    for(nclass = 0; nclass < cm->numClasses; nclass++) {
      if (grid->dim == 1) {
        ierr = PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_1D_CONSTANT,  &isconst);    CHKERRQ(ierr);
        ierr = PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_1D_LINEAR,    &islin);      CHKERRQ(ierr);
        ierr = PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_1D_QUADRATIC, &isquad);     CHKERRQ(ierr);
      } else if (grid->dim == 2) {
        ierr = PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_2D_LINEAR,    &islin);      CHKERRQ(ierr);
        ierr = PetscTypeCompare((PetscObject) disc, DISCRETIZATION_TRIANGULAR_2D_QUADRATIC, &isquad);     CHKERRQ(ierr);
      }
      /* These only work for pairs in 2D (all are fine in 1D, how do I generalize this?) */
      if (isconst == PETSC_TRUE) {
        cm->fieldClasses[f][nclass] = (nclass == 1 ? 1 : 0);
      } else if (islin == PETSC_TRUE) {
        cm->fieldClasses[f][nclass] = (nclass == 0 ? 1 : 0);
      } else if (isquad == PETSC_TRUE) {
        cm->fieldClasses[f][nclass] = 1;
      } else {
        SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid discretization type");
      }
    }
  }

  /* Setup classes */
  ierr = PetscMemzero(cm->classes, cm->numOverlapNodes * sizeof(int));                                    CHKERRQ(ierr);
  if (cm->numClasses > 1) {
    /* Set all midnodes to class 1 */
    ierr = PartitionGetNumOverlapElements(part, &numOverlapElements);                                     CHKERRQ(ierr);
    ierr = MeshGetNumCorners(mesh, &numCorners);                                                          CHKERRQ(ierr);
    for(elem = 0; elem < numOverlapElements; elem++) {
      for(corner = grid->dim+1; corner < numCorners; corner++) {
        ierr = MeshGetNodeFromElement(mesh, elem, corner, &node);                                         CHKERRQ(ierr);
        if (cm->classes[node] == 0)
          cm->classes[node] = 1;
      }
    }
  }

  /* Communicate ghost classes */
  ierr = PartitionGhostNodeExchange(part, INSERT_VALUES, SCATTER_FORWARD, cm->classes, &cm->classes[cm->numNodes]);
  CHKERRQ(ierr);

  /* Setup class sizes */
  ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizes);                                      CHKERRQ(ierr);
  ierr = PetscMemzero(cm->classSizes, cm->numClasses * sizeof(int));                                      CHKERRQ(ierr);
  for(nclass = 0; nclass < cm->numClasses; nclass++) {
    for(f = 0; f < cm->numFields; f++) {
      field = cm->fields[f];
      if (cm->fieldClasses[f][nclass]) {
        ierr = GridGetFieldComponents(grid, field, &comp);                                                CHKERRQ(ierr);
        cm->classSizes[nclass] += comp;
      }
    }
  }

  /* Initialize constraint support */
  ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->isClassConstrained);                              CHKERRQ(ierr);
  ierr = PetscMalloc(cm->numClasses * sizeof(int), &cm->classSizeDiffs);                                  CHKERRQ(ierr);
  ierr = PetscMemzero(cm->isClassConstrained, cm->numClasses * sizeof(int));                              CHKERRQ(ierr);
  ierr = PetscMemzero(cm->classSizeDiffs,     cm->numClasses * sizeof(int));                              CHKERRQ(ierr);

  PetscLogObjectMemory(cm, (cm->numFields + cm->numOverlapNodes + (cm->numFields + 4)*cm->numClasses) * sizeof(int) +
                   cm->numFields * sizeof(int *));
  PetscFunctionReturn(0);
}
EXTERN_C_END