#define PETSCDM_DLL
#include <petsc-private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/

#undef __FUNCT__
#define __FUNCT__ "DMPlexCreateFluentFromFile"
/*@C
  DMPlexCreateFluentFromFile - Create a DMPlex mesh from a Fluent mesh file

+ comm        - The MPI communicator
. filename    - Name of the Fluent mesh file
- interpolate - Create faces and edges in the mesh

  Output Parameter:
. dm  - The DM object representing the mesh

  Level: beginner

.seealso: DMPlexCreateFromFile(), DMPlexCreateFluent(), DMPlexCreate()
@*/
PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
{
  PetscViewer     viewer;
  PetscErrorCode  ierr;

  PetscFunctionBegin;
  /* Create file viewer and build plex */
  ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr);
  ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr);
  ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr);
  ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr);
  ierr = DMPlexCreateFluent(comm, viewer, interpolate, dm);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DMPlexCreateFluent_ReadString"
PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
{
  PetscInt i = 0;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  do {ierr = PetscViewerRead(viewer, &(buffer[i++]), 1, PETSC_CHAR);CHKERRQ(ierr);}
  while (buffer[i-1] != '\0' && buffer[i-1] != delim);
  buffer[i] = '\0';
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DMPlexCreateFluent_ReadSection"
PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
{
  char           buffer[PETSC_MAX_PATH_LEN];
  int            snum;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  /* Fast-forward to next section and derive its index */
  ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
  ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ' ');CHKERRQ(ierr);
  snum = sscanf(buffer, "%d", &(s->index));
  /* If we can't match an index return -1 to signal end-of-file */
  if (snum < 1) {s->index = -1;   PetscFunctionReturn(0);}

  if (s->index == 0) {           /* Comment */
    ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);

  } else if (s->index == 2) {    /* Dimension */
    ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
    snum = sscanf(buffer, "%d", &(s->nd));
    if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");

  } else if (s->index == 10) {   /* Vertices */
    ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
    snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
    if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    if (s->zoneID > 0) {
      PetscInt numCoords = s->last - s->first + 1;
      ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
      ierr = PetscMalloc1(s->nd*numCoords, (PetscScalar**)&s->data);CHKERRQ(ierr);
      ierr = PetscViewerRead(viewer, (PetscScalar*)s->data, numCoords*s->nd, PETSC_REAL);CHKERRQ(ierr);
    }

  } else if (s->index == 12) {   /* Cells */
    ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
    snum = sscanf(buffer, "(%x", &(s->zoneID));
    if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    if (s->zoneID == 0) {  /* Header section */
      snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
      if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    } else {               /* Data section */
      snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
      if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    }
    ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);

  } else if (s->index == 13) {   /* Faces */
    ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
    snum = sscanf(buffer, "(%x", &(s->zoneID));
    if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    if (s->zoneID == 0) {  /* Header section */
      snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
      if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    } else {               /* Data section */
      PetscInt f, e, numEntries, numFaces;
      snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
      if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
      ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
      switch (s->nd) {
      case 0: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed faces in Fluent files are not supported");
      case 2: numEntries = 2 + 2; break;  /* linear */
      case 3: numEntries = 2 + 3; break;  /* triangular */
      case 4: numEntries = 2 + 4; break;  /* quadrilateral */
      default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
      }
      numFaces = s->last-s->first + 1;
      ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr);
      for (f = 0; f < numFaces; f++) {
        for (e = 0; e < numEntries; e++) {
          ierr = PetscViewerRead(viewer, buffer, 1, PETSC_STRING);CHKERRQ(ierr);
          snum = sscanf(buffer, "%x", &(((PetscInt*)s->data)[f*numEntries + e]));
          if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
        }
      }
      ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
    }
    ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);

  } else {                       /* Unknown section type */
    PetscInt depth = 1;
    do {
      /* Match parentheses when parsing unknown sections */
      do {ierr = PetscViewerRead(viewer, &(buffer[0]), 1, PETSC_CHAR);CHKERRQ(ierr);}
      while (buffer[0] != '(' && buffer[0] != ')');
      if (buffer[0] == '(') depth++;
      if (buffer[0] == ')') depth--;
    } while (depth > 0);
    ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '\n');CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DMPlexCreateFluent"
/*@C
  DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file.

  Collective on comm

  Input Parameters:
+ comm  - The MPI communicator
. viewer - The Viewer associated with a Fluent mesh file
- interpolate - Create faces and edges in the mesh

  Output Parameter:
. dm  - The DM object representing the mesh

  Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm

  Level: beginner

.keywords: mesh, fluent, case
.seealso: DMPLEX, DMCreate()
@*/
PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
{
  PetscMPIInt    rank;
  PetscInt       c, f, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
  PetscInt       numFaces = PETSC_DETERMINE, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
  PetscInt      *faces = NULL, *cellVertices;
  PetscInt       d, coordSize;
  PetscScalar   *coords, *coordsIn = NULL;
  PetscSection   coordSection;
  Vec            coordinates;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);

  if (!rank) {
    FluentSection s;
    do {
      ierr = DMPlexCreateFluent_ReadSection(viewer, &s);CHKERRQ(ierr);
      if (s.index == 2) {                 /* Dimension */
        dim = s.nd;

      } else if (s.index == 10) {         /* Vertices */
        if (s.zoneID == 0) numVertices = s.last;
        else {
          if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
          coordsIn = s.data;
        }

      } else if (s.index == 12) {         /* Cells */
        if (s.zoneID == 0) numCells = s.last;
        else {
          switch (s.nd) {
          case 0: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed elements in Fluent files are not supported");
          case 1: numCellVertices = 3; break;  /* triangular */
          case 2: numCellVertices = 4; break;  /* tetrahedral */
          case 3: numCellVertices = 4; break;  /* quadrilateral */
          case 4: numCellVertices = 8; break;  /* hexahedral */
          case 5: numCellVertices = 5; break;  /* pyramid */
          case 6: numCellVertices = 6; break;  /* wedge */
          default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell element-type in Fluent file");
          }
        }

      } else if (s.index == 13) {         /* Facets */
        if (s.zoneID == 0) {  /* Header section */
          numFaces = s.last - s.first + 1;
        } else {              /* Data section */
          if (s.nd == 0 || (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices)) {
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are currently not supported");
          }
          if (numFaces == PETSC_DETERMINE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
          if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
          numFaceEntries = numFaceVertices + 2;
          if (!faces) {ierr = PetscMalloc1(numFaces*numFaceEntries, &faces);CHKERRQ(ierr);}
          ierr = PetscMemcpy(&(faces[(s.first-1)*numFaceEntries]), s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));CHKERRQ(ierr);
          ierr = PetscFree(s.data);CHKERRQ(ierr);
        }
      }
    } while (s.index >= 0);
  }
  ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
  if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");

  /* Allocate cell-vertex mesh */
  ierr = DMCreate(comm, dm);CHKERRQ(ierr);
  ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
  ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
  ierr = DMPlexSetChart(*dm, 0, numCells + numVertices);CHKERRQ(ierr);
  if (!rank) {
    if (numCellVertices == PETSC_DETERMINE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type not specified in Fluent file");
    for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(*dm, c, numCellVertices);CHKERRQ(ierr);}
  }
  ierr = DMSetUp(*dm);CHKERRQ(ierr);

  if (!rank) {
    /* Derive cell-vertex list from face-vertex and face-cell maps */
    if (numCells == PETSC_DETERMINE || numCellVertices == PETSC_DETERMINE) {
      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Insufficent cell header information in Fluent file");
    }
    ierr = PetscMalloc1(numCells*numCellVertices, &cellVertices);CHKERRQ(ierr);
    for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1;
    for (f = 0; f < numFaces; f++) {
      PetscInt *cell;
      const PetscInt cl = faces[f*numFaceEntries + numFaceVertices];
      const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1];
      const PetscInt *face = &(faces[f*numFaceEntries]);

      if (cl > 0) {
        cell = &(cellVertices[(cl-1) * numCellVertices]);
        for (v = 0; v < numFaceVertices; v++) {
          PetscBool found = PETSC_FALSE;
          for (c = 0; c < numCellVertices; c++) {
            if (cell[c] < 0) break;
            if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
          }
          if (!found) cell[c] = face[v]-1 + numCells;
        }
      }
      if (cr > 0) {
        cell = &(cellVertices[(cr-1) * numCellVertices]);
        for (v = 0; v < numFaceVertices; v++) {
          PetscBool found = PETSC_FALSE;
          for (c = 0; c < numCellVertices; c++) {
            if (cell[c] < 0) break;
            if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
          }
          if (!found) cell[c] = face[v]-1 + numCells;
        }
      }
    }
    for (c = 0; c < numCells; c++) {
      ierr = DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));CHKERRQ(ierr);
    }
  }
  ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
  ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
  if (interpolate) {
    DM idm = NULL;

    ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
    ierr = DMDestroy(dm);CHKERRQ(ierr);
    *dm  = idm;
  }

  /* Read coordinates */
  ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
  ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
  ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
  ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr);
  for (v = numCells; v < numCells+numVertices; ++v) {
    ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
    ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
  }
  ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
  ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
  ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
  ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
  ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr);
  ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
  if (!rank) {
    for (v = 0; v < numVertices; ++v) {
      for (d = 0; d < dim; ++d) {
        coords[v*dim+d] = coordsIn[v*dim+d];
      }
    }
  }
  ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
  ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
  ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
  if (!rank) {
    ierr = PetscFree(cellVertices);CHKERRQ(ierr);
    ierr = PetscFree(faces);CHKERRQ(ierr);
    ierr = PetscFree(coordsIn);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}