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

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

#undef  __FUNCT__
#define __FUNCT__ "GridCreate"
/*@ 
  GridCreate - This function creates an empty grid. The type can then be set with GridSetType().

  Collective on Mesh

  Input Parameter:
. mesh - The Mesh

  Output Parameter:
. grid - The grid

  Options Database Keys:
. -grid_explicit_constraints - Constraints are implemented at the element level instead of through projectors
. -grid_int_type <type>      - The interpolation type, e.g local, l2, etc.

  Level: beginner

.keywords: grid, create
.seealso: GridSetType(), GridSetUp(), GridDestroy(), MeshCreate()
@*/
int GridCreate(Mesh mesh, Grid *grid) {
  MPI_Comm comm;
  Grid     g;
  int      dim;
  int      ierr;

  PetscFunctionBegin;
  PetscValidPointer(grid);
  *grid = PETSC_NULL;
#ifndef PETSC_USE_DYNAMIC_LIBRARIES
  ierr = GridInitializePackage(PETSC_NULL);                                                               CHKERRQ(ierr);
#endif

  ierr = PetscObjectGetComm((PetscObject) mesh, &comm);                                                   CHKERRQ(ierr);
  PetscHeaderCreate(g, _Grid, struct _GridOps, GRID_COOKIE, -1, "Grid", comm, GridDestroy, GridView);
  PetscLogObjectCreate(g);
  PetscLogObjectMemory(g, sizeof(struct _Grid));
  ierr = PetscMemzero(g->ops, sizeof(struct _GridOps));                                                   CHKERRQ(ierr);
  g->bops->publish    = PETSC_NULL /* GridPublish_Petsc */;
  g->type_name        = PETSC_NULL;
  g->serialize_name   = PETSC_NULL;

  /* General grid description */
  ierr = MeshGetDimension(mesh, &dim);                                                                    CHKERRQ(ierr);
  g->dim           = dim;
  g->mesh          = mesh;
  g->gridparent    = PETSC_NULL;
  g->setupcalled   = PETSC_FALSE;
  g->bdSetupCalled = PETSC_FALSE;
  g->data          = PETSC_NULL;
  g->usr           = PETSC_NULL;
  PetscObjectReference((PetscObject) mesh);

  /* Field variables */
  g->numFields  = 0;
  g->maxFields  = 1;
  ierr = PetscMalloc(g->maxFields * sizeof(Field), &g->fields);                                           CHKERRQ(ierr);

  /* Class structure */
  g->cm = PETSC_NULL;

  /* Default variable orderings */
  g->order    = PETSC_NULL;
  g->locOrder = PETSC_NULL;

  /* Ghost variable scatter */
  g->ghostVec     = PETSC_NULL;
  g->ghostScatter = PETSC_NULL;

  /* Constraint variables */
  g->isConstrained       = PETSC_FALSE;
  g->explicitConstraints = PETSC_FALSE;
  g->numNewFields        = -1;
  g->constraintCM        = PETSC_NULL;
  g->constraintOrder     = PETSC_NULL;
  g->constraintOrdering  = PETSC_NULL;
  g->constraintMatrix    = PETSC_NULL;
  g->constraintCtx       = PETSC_NULL;

  /* Problem variables */
  g->numRhsFuncs = 0;
  g->maxRhsFuncs = 1;
  ierr = PetscMalloc(g->maxRhsFuncs * sizeof(GridFunc), &g->rhsFuncs);                                    CHKERRQ(ierr);
  g->numRhsOps   = 0;
  g->maxRhsOps   = 1;
  ierr = PetscMalloc(g->maxRhsOps   * sizeof(GridOp),   &g->rhsOps);                                      CHKERRQ(ierr);
  g->numMatOps   = 0;
  g->maxMatOps   = 1;
  ierr = PetscMalloc(g->maxMatOps   * sizeof(GridOp),   &g->matOps);                                      CHKERRQ(ierr);

  /* Problem query variables */
  g->activeMatOp       = -1;
  g->activeNonlinearOp = -1;

  /* Assembly variables */
  g->numActiveFields  = 0;
  g->maxActiveFields  = 1;
  ierr = PetscMalloc(g->maxActiveFields * sizeof(int), &g->defaultFields);                                CHKERRQ(ierr);
  g->vec              = PETSC_NULL;
  g->mat              = PETSC_NULL;
  g->ghostElementVec  = PETSC_NULL;
  g->ALEActive        = PETSC_FALSE;
  g->activeOpTypes[0] = PETSC_TRUE;
  g->activeOpTypes[1] = PETSC_TRUE;
  g->activeOpTypes[2] = PETSC_TRUE;

  /* Boundary condition variables */
  g->reduceSystem      = PETSC_FALSE;
  g->reduceElement     = PETSC_FALSE;
  g->reduceElementArgs = PETSC_FALSE;
  g->reductionCM       = PETSC_NULL;
  g->reduceOrder       = PETSC_NULL;
  g->locReduceOrder    = PETSC_NULL;
  g->bdReduceVec       = PETSC_NULL;
  g->bdReduceVecOld    = PETSC_NULL;
  g->bdReduceVecDiff   = PETSC_NULL;
  g->bdReduceVecCur    = PETSC_NULL;
  g->bdReduceMat       = PETSC_NULL;
  g->reduceVec         = PETSC_NULL;
  g->reduceAlpha       = 1.0;
  g->reduceContext     = PETSC_NULL;
  g->numBC             = 0;
  g->maxBC             = 1;
  ierr = PetscMalloc(g->maxBC      * sizeof(GridBC), &g->bc);                                             CHKERRQ(ierr);
  g->numPointBC        = 0;
  g->maxPointBC        = 1;
  ierr = PetscMalloc(g->maxPointBC * sizeof(GridBC), &g->pointBC);                                        CHKERRQ(ierr);

  /* Boundary iteration variables */
  ierr = MeshGetNumBoundaries(mesh, &g->numBd);                                                           CHKERRQ(ierr);
  ierr = PetscMalloc(g->numBd * sizeof(int *), &g->bdSize);                                               CHKERRQ(ierr);
  PetscLogObjectMemory(g, g->numBd * sizeof(int *));
  ierr = PetscMemzero(g->bdSize, g->numBd * sizeof(int *));                                               CHKERRQ(ierr);
  g->bdOrder    = PETSC_NULL;
  g->bdLocOrder = PETSC_NULL;

  /* Setup matrix-free */
  g->isMatrixFree      = PETSC_FALSE;
  g->matrixFreeArg     = PETSC_NULL;
  g->matrixFreeContext = PETSC_NULL;

  /* Interpolation variables */
  g->interpolationType = INTERPOLATION_LOCAL;

  /* Graphics extras */
  g->viewField = -1;
  g->viewComp  = 0;

  *grid = g;
  PetscFunctionReturn(0);
}

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

  Collective on MPI_Comm

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

  Output Parameter:
. grid   - The grid

  Level: beginner

.keywords: grid, serialize
.seealso: PartitionSerialize(), GridSerialize()
@*/
int GridSerialize(MPI_Comm comm, Grid *grid, PetscViewer viewer, PetscTruth store)
{
  int      (*serialize)(MPI_Comm, Grid *, PetscViewer, PetscTruth);
  int        fd, len;
  char      *name;
  PetscTruth match;
  int        ierr;

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

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

  if (!GridSerializeRegisterAllCalled) {
    ierr = GridSerializeRegisterAll(PETSC_NULL);                                                          CHKERRQ(ierr);
  }
  if (!GridSerializeList) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Could not find table of methods");

  if (store) {
    PetscValidHeaderSpecific(*grid, GRID_COOKIE);
    ierr = PetscStrlen((*grid)->class_name, &len);                                                        CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &len,                     1,   PETSC_INT,  0);                            CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  (*grid)->class_name,     len, PETSC_CHAR, 0);                            CHKERRQ(ierr);
    ierr = PetscStrlen((*grid)->serialize_name, &len);                                                    CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &len,                     1,   PETSC_INT,  0);                            CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd,  (*grid)->serialize_name, len, PETSC_CHAR, 0);                            CHKERRQ(ierr);
    ierr = PetscFListFind(comm, GridSerializeList, (*grid)->serialize_name, (void (**)(void)) &serialize); CHKERRQ(ierr);
    if (!serialize) SETERRQ(PETSC_ERR_ARG_WRONG, "Type cannot be serialized");
    ierr = (*serialize)(comm, grid, viewer, store);                                                       CHKERRQ(ierr);
  } else {
    ierr = PetscBinaryRead(fd, &len,    1,   PETSC_INT);                                                  CHKERRQ(ierr);
    ierr = PetscMalloc((len+1) * sizeof(char), &name);                                                    CHKERRQ(ierr);
    name[len] = 0;
    ierr = PetscBinaryRead(fd,  name,   len, PETSC_CHAR);                                                 CHKERRQ(ierr);
    ierr = PetscStrcmp(name, "Grid", &match);                                                             CHKERRQ(ierr);
    ierr = PetscFree(name);                                                                               CHKERRQ(ierr);
    if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONG, "Non-grid object");
    /* Dispatch to the correct routine */
    ierr = PetscBinaryRead(fd, &len,    1,   PETSC_INT);                                                  CHKERRQ(ierr);
    ierr = PetscMalloc((len+1) * sizeof(char), &name);                                                    CHKERRQ(ierr);
    name[len] = 0;
    ierr = PetscBinaryRead(fd,  name,   len, PETSC_CHAR);                                                 CHKERRQ(ierr);
    ierr = PetscFListFind(comm, GridSerializeList, name, (void (**)(void)) &serialize);                   CHKERRQ(ierr);
    if (!serialize) SETERRQ(PETSC_ERR_ARG_WRONG, "Type cannot be serialized");
    ierr = (*serialize)(comm, grid, viewer, store);                                                       CHKERRQ(ierr);
    ierr = PetscStrfree((*grid)->serialize_name);                                                         CHKERRQ(ierr);
    (*grid)->serialize_name = name;
  }

  PetscFunctionReturn(0);
}