static char help[] = "Test for interpolating a projected field.\n\n\n";

#include <petscdmplex.h>
#include <petscds.h>

typedef struct {
  /* Domain and mesh definition */
  PetscBool interpolate;         /* Generate intermediate mesh elements, e.g. faces and edges */
  PetscBool reversecell;         /* Use DMPlexReverseCell */ 
  PetscInt  maxprojectionheight; /* The max projection height to use */
} AppCtx;

PetscErrorCode coordinates_function(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
  const PetscInt Ncomp = dim;

  PetscInt       comp;
  for (comp = 0; comp < Ncomp; ++comp) u[comp] = x[comp];
  return 0;
}

void coordinates_field(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    PetscReal t, const PetscReal x[], PetscScalar f0[])
{
  const PetscInt Ncomp = dim;

  PetscInt       comp;
  for (comp = 0; comp < Ncomp; ++comp) f0[comp] = u[comp];
}


#undef __FUNCT__
#define __FUNCT__ "ProcessOptions"
PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
  options->interpolate         = PETSC_FALSE;
  options->reversecell         = PETSC_FALSE;
  options->maxprojectionheight = 0;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = PetscOptionsBegin(comm, "", "Reverse cell test options", "DMPLEX");CHKERRQ(ierr);
  ierr = PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "testreverse.c", options->interpolate, &options->interpolate, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-reversecell", "Use DMPlexReverseCell", "testreverse.c", options->reversecell, &options->reversecell, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsInt("-maxprojectionheight", "Maximum projection height", "testreverse.c", options->maxprojectionheight, &options->maxprojectionheight, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CreateMesh"
PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
{
  PetscInt       dim                   = 3;
  PetscBool      interpolate           = user->interpolate, reversecell = user->reversecell;
  const int      cells[]               = {0, 1, 2, 3};
  const double   reverseVertexCoords[] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0};
  const double   vertexCoords[]        = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0};
  PetscInt       numCells = 1, numVertices = 4, numCorners = 4;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, reversecell ? reverseVertexCoords : vertexCoords, dm);CHKERRQ(ierr);
  {
    PetscReal v0[3], J[9], invJ[9], detJ;

    ierr = DMPlexComputeCellGeometryFEM(*dm, 0, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
    if (detJ < 0) {
      ierr = DMPlexReverseCell(*dm, 0);CHKERRQ(ierr);
      ierr = PetscPrintf(comm, "Cell reversed\n");CHKERRQ(ierr);
    }
  }
  ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) *dm, "Simplicial Mesh");CHKERRQ(ierr);
  ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "SetupDiscretization"
PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
{
  const PetscInt  dim = 3;
  PetscInt        f;
  const PetscBool simplex = PETSC_TRUE;
  PetscFE         fe[3];
  PetscDS         prob, probAux;
  PetscErrorCode  ierr;

  PetscFunctionBeginUser;
  /* Create finite element */
  ierr = PetscFECreateDefault(dm, dim, dim, simplex, "projectedcoordinates_", -1, &fe[0]);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe[0], "projectedcoordinates");CHKERRQ(ierr);
  ierr = PetscFECreateDefault(dm, dim, dim, simplex, "interpolatedcoordinates_", -1, &fe[1]);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe[1], "interpolatedcoordinates");CHKERRQ(ierr);
  ierr = PetscFECreateDefault(dm, dim, dim, simplex, "interpolatedprojectedcoordinates_", -1, &fe[2]);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe[2], "interpolatedprojectedcoordinates");CHKERRQ(ierr);

  for (f = 1; f < 3; f++) {
    PetscReal      *points, *weights;
    PetscDualSpace  sp;
    PetscQuadrature q;

    ierr = PetscFEGetDualSpace(fe[1], &sp);CHKERRQ(ierr);
    ierr = PetscDualSpaceSetType(sp, PETSCDUALSPACESIMPLE);CHKERRQ(ierr);
    ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &q);CHKERRQ(ierr);
    ierr = PetscMalloc1(dim, &points);CHKERRQ(ierr);
    ierr = PetscMalloc1(1, &weights);CHKERRQ(ierr);
    weights[0] = 1.0;
    ierr = DMPlexComputeCellGeometryFVM(dm, 0, NULL, points, NULL);CHKERRQ(ierr);
    ierr = PetscQuadratureSetData(q, dim, 1, points, weights);CHKERRQ(ierr);
    ierr = PetscDualSpaceSimpleSetDimension(sp, 1);CHKERRQ(ierr);
    ierr = PetscDualSpaceSimpleSetFunctional(sp, 0, q);CHKERRQ(ierr);
  }


  ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
  ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe[0]);CHKERRQ(ierr);
  ierr = PetscDSSetDiscretization(prob, 1, (PetscObject) fe[1]);CHKERRQ(ierr);
  ierr = PetscDSSetDiscretization(prob, 2, (PetscObject) fe[2]);CHKERRQ(ierr);

  ierr = PetscDSViewFromOptions(prob, NULL, "-ds_view");CHKERRQ(ierr);

  ierr = PetscFEDestroy(&fe[0]);CHKERRQ(ierr);
  ierr = PetscFEDestroy(&fe[1]);CHKERRQ(ierr);
  ierr = PetscFEDestroy(&fe[2]);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv)
{
  PetscInt       numFields, f;
  char         **fieldNames;
  IS            *fields;
  DM             dm;                   
  Vec            u, v, subvec;
  AppCtx         user;                 
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
  ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
  ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr);

  ierr = SetupDiscretization(dm, &user);CHKERRQ(ierr);
  ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr);

  ierr = DMPlexSetMaxProjectionHeight(dm, user.maxprojectionheight);CHKERRQ(ierr);

  ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr);

  ierr = DMCreateFieldDecomposition(dm, &numFields, &fieldNames, &fields, NULL);CHKERRQ(ierr);

  /* Project the coordinates on the projected coordinate field and on the interpolated coordinate field */
  PetscErrorCode (*funcs_function[])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void* ctx) = {coordinates_function, coordinates_function, NULL};
  void *ctxs[] = {&user, &user, &user};
  ierr = DMProjectFunction(dm, 0.0, funcs_function, ctxs, INSERT_VALUES, u);CHKERRQ(ierr);

  /* View the projected coordinates field */
  ierr = VecGetSubVector(u, fields[0], &subvec);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) subvec, fieldNames[0]);CHKERRQ(ierr);
  ierr = VecViewFromOptions(subvec, NULL, "-projectedcoordinates_vec_view");CHKERRQ(ierr);
  ierr = VecRestoreSubVector(u, fields[0], &subvec);CHKERRQ(ierr);

  /* View the interpolated coordinates field */
  ierr = VecGetSubVector(u, fields[1], &subvec);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) subvec, fieldNames[1]);CHKERRQ(ierr);
  ierr = VecViewFromOptions(subvec, NULL, "-interpolatedcoordinates_vec_view");CHKERRQ(ierr);
  ierr = VecRestoreSubVector(u, fields[1], &subvec);CHKERRQ(ierr);

  ierr = DMPlexSetMaxProjectionHeight(dm, 0);CHKERRQ(ierr);
  ierr = VecDuplicate(u, &v);CHKERRQ(ierr);
  /* Use the projected coordinate field to interpolate to the interpolated projected coordinate field */
  void (*funcs_field[])(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscScalar[]) = {NULL, NULL, coordinates_field};
  ierr = DMProjectField(dm, u, funcs_field, INSERT_VALUES, v);CHKERRQ(ierr);
  
  /* View the interpolated projected coordinates field */
  ierr = VecGetSubVector(v, fields[2], &subvec);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) subvec, fieldNames[2]);CHKERRQ(ierr);
  ierr = VecViewFromOptions(subvec, NULL, "-interpolatedprojectedcoordinates_vec_view");CHKERRQ(ierr);
  ierr = VecRestoreSubVector(v, fields[2], &subvec);CHKERRQ(ierr);

  for (f = 0; f < numFields; ++f) {
    ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr);
    ierr = ISDestroy(&fields[f]);CHKERRQ(ierr);
  }
  ierr = PetscFree(fieldNames);CHKERRQ(ierr);
  ierr = PetscFree(fields);CHKERRQ(ierr);

  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&v);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}