static char help[] = "Tests DMDA ghost coordinates\n\n";

#include <petscdm.h>
#include <petscdmda.h>

static PetscErrorCode CompareGhostedCoords(Vec gc1, Vec gc2)
{
  PetscReal nrm, gnrm;
  Vec       tmp;

  PetscFunctionBeginUser;
  PetscCall(VecDuplicate(gc1, &tmp));
  PetscCall(VecWAXPY(tmp, -1.0, gc1, gc2));
  PetscCall(VecNorm(tmp, NORM_INFINITY, &nrm));
  PetscCallMPI(MPIU_Allreduce(&nrm, &gnrm, 1, MPIU_REAL, MPIU_MAX, PETSC_COMM_WORLD));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "norm of difference of ghosted coordinates %8.2e\n", (double)gnrm));
  PetscCall(VecDestroy(&tmp));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TestQ2Q1DA(void)
{
  DM       Q2_da, Q1_da, cda;
  PetscInt mx, my, mz;
  Vec      coords, gcoords, gcoords2;

  PetscFunctionBeginUser;
  mx = 7;
  my = 11;
  mz = 13;
  PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, mz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 3, 2, 0, 0, 0, &Q2_da));
  PetscCall(DMSetFromOptions(Q2_da));
  PetscCall(DMSetUp(Q2_da));
  PetscCall(DMDASetUniformCoordinates(Q2_da, -1.0, 1.0, -2.0, 2.0, -3.0, 3.0));
  PetscCall(DMGetCoordinates(Q2_da, &coords));
  PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, mz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 3, 1, 0, 0, 0, &Q1_da));
  PetscCall(DMSetFromOptions(Q1_da));
  PetscCall(DMSetUp(Q1_da));
  PetscCall(DMSetCoordinates(Q1_da, coords));

  /* Get ghost coordinates one way */
  PetscCall(DMGetCoordinatesLocal(Q1_da, &gcoords));

  /* And another */
  PetscCall(DMGetCoordinates(Q1_da, &coords));
  PetscCall(DMGetCoordinateDM(Q1_da, &cda));
  PetscCall(DMGetLocalVector(cda, &gcoords2));
  PetscCall(DMGlobalToLocalBegin(cda, coords, INSERT_VALUES, gcoords2));
  PetscCall(DMGlobalToLocalEnd(cda, coords, INSERT_VALUES, gcoords2));

  PetscCall(CompareGhostedCoords(gcoords, gcoords2));
  PetscCall(DMRestoreLocalVector(cda, &gcoords2));

  PetscCall(VecScale(coords, 10.0));
  PetscCall(VecScale(gcoords, 10.0));
  PetscCall(DMGetCoordinatesLocal(Q1_da, &gcoords2));
  PetscCall(CompareGhostedCoords(gcoords, gcoords2));

  PetscCall(DMDestroy(&Q2_da));
  PetscCall(DMDestroy(&Q1_da));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char **argv)
{
  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, 0, help));
  PetscCall(TestQ2Q1DA());
  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

   test:
      nsize: 2

TEST*/