#include <petsc.h>

/* Run the test for 2D using

     -dm_plex_simplex 0 -dm_plex_box_faces 16,16 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 -dm_plex_box_bd periodic,periodic \
       -dm_plex_periodic_cut -dm_view hdf5:sol.h5 -vec_view hdf5:sol.h5::append -dm_plex_create_fv_ghost_cells

   and for 3D,

     -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 16,16,16 -dm_plex_box_lower -5,-5,-5 -dm_plex_box_upper 5,5,5 \
       -dm_plex_box_bd periodic,periodic,periodic -dm_plex_periodic_cut -dm_view hdf5:sol.h5 -vec_view hdf5:sol.h5::append
*/

static PetscErrorCode vortex(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
{
  const PetscReal gamma = 1.4;
  const PetscReal beta  = 5.0;
  const PetscReal rsq   = PetscSqr(x[0]) + PetscSqr(x[1]);

  PetscFunctionBeginUser;
  /* Set density */
  u[0] = PetscPowReal(1. - ((gamma - 1.)*beta*beta)/(8.0*gamma*PETSC_PI*PETSC_PI) * PetscExpReal(1. - rsq), 1./(gamma - 1.));
  PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
  DM             dm;
  PetscDS        ds;
  PetscFV        fvm;
  Vec            sol;
  PetscSimplePointFunc func[1] = {vortex};
  PetscInt       dim, Nf = 1;
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv, NULL, NULL);if (ierr) return ierr;
  ierr = DMCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr);
  ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
  ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
  ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
  ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
  /* TODO Add ghost cells */

  ierr = PetscFVCreate(PETSC_COMM_WORLD, &fvm);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fvm, "FV solver");CHKERRQ(ierr);
  ierr = PetscFVSetFromOptions(fvm);CHKERRQ(ierr);
  ierr = PetscFVSetNumComponents(fvm, Nf);CHKERRQ(ierr);
  ierr = PetscFVSetSpatialDimension(fvm, dim);CHKERRQ(ierr);
  ierr = PetscFVSetComponentName(fvm, 0, "Density");CHKERRQ(ierr);

  ierr = DMAddField(dm, NULL, (PetscObject) fvm);CHKERRQ(ierr);
  ierr = PetscFVDestroy(&fvm);CHKERRQ(ierr);
  ierr = DMCreateDS(dm);CHKERRQ(ierr);
  ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);
  ierr = PetscDSSetFromOptions(ds);CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(dm, &sol);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) sol, "Solution");CHKERRQ(ierr);
  ierr = VecZeroEntries(sol);CHKERRQ(ierr);
  ierr = DMProjectFunction(dm, 0.0, func, NULL, INSERT_ALL_VALUES, sol);CHKERRQ(ierr);
  ierr = VecViewFromOptions(sol, NULL, "-vec_view");CHKERRQ(ierr);

  ierr = VecDestroy(&sol);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}