static char help[] = "Checks the functionality of DMGetInterpolation() on deformed grids.\n\n";

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

typedef struct _n_CCmplx CCmplx;
struct _n_CCmplx {
  PetscReal real;
  PetscReal imag;
};

CCmplx CCmplxPow(CCmplx a, PetscReal n)
{
  CCmplx    b;
  PetscReal r, theta;
  r      = PetscSqrtReal(a.real * a.real + a.imag * a.imag);
  theta  = PetscAtan2Real(a.imag, a.real);
  b.real = PetscPowReal(r, n) * PetscCosReal(n * theta);
  b.imag = PetscPowReal(r, n) * PetscSinReal(n * theta);
  return b;
}
CCmplx CCmplxExp(CCmplx a)
{
  CCmplx b;
  b.real = PetscExpReal(a.real) * PetscCosReal(a.imag);
  b.imag = PetscExpReal(a.real) * PetscSinReal(a.imag);
  return b;
}
CCmplx CCmplxSqrt(CCmplx a)
{
  CCmplx    b;
  PetscReal r, theta;
  r      = PetscSqrtReal(a.real * a.real + a.imag * a.imag);
  theta  = PetscAtan2Real(a.imag, a.real);
  b.real = PetscSqrtReal(r) * PetscCosReal(0.5 * theta);
  b.imag = PetscSqrtReal(r) * PetscSinReal(0.5 * theta);
  return b;
}
CCmplx CCmplxAdd(CCmplx a, CCmplx c)
{
  CCmplx b;
  b.real = a.real + c.real;
  b.imag = a.imag + c.imag;
  return b;
}
PetscScalar CCmplxRe(CCmplx a)
{
  return a.real;
}
PetscScalar CCmplxIm(CCmplx a)
{
  return a.imag;
}

PetscErrorCode DAApplyConformalMapping(DM da, PetscInt idx)
{
  PetscInt     i, n;
  PetscInt     sx, nx, sy, ny, sz, nz, dim;
  Vec          Gcoords;
  PetscScalar *XX;
  PetscScalar  xx, yy, zz;
  DM           cda;

  PetscFunctionBeginUser;
  if (idx == 0) {
    PetscFunctionReturn(PETSC_SUCCESS);
  } else if (idx == 1) { /* dam break */
    PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
  } else if (idx == 2) { /* stagnation in a corner */
    PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, 0.0, 1.0, -1.0, 1.0));
  } else if (idx == 3) { /* nautilis */
    PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
  } else if (idx == 4) PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));

  PetscCall(DMGetCoordinateDM(da, &cda));
  PetscCall(DMGetCoordinates(da, &Gcoords));

  PetscCall(VecGetArray(Gcoords, &XX));
  PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
  PetscCall(DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
  PetscCall(VecGetLocalSize(Gcoords, &n));
  n = n / dim;

  for (i = 0; i < n; i++) {
    if ((dim == 3) && (idx != 2)) {
      PetscScalar Ni[8];
      PetscScalar xi   = XX[dim * i];
      PetscScalar eta  = XX[dim * i + 1];
      PetscScalar zeta = XX[dim * i + 2];
      PetscScalar xn[] = {-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0};
      PetscScalar yn[] = {-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0};
      PetscScalar zn[] = {-0.1, -4.0, -0.2, -1.0, 0.1, 4.0, 0.2, 1.0};
      PetscInt    p;

      Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
      Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
      Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
      Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);

      Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
      Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
      Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
      Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);

      xx = yy = zz = 0.0;
      for (p = 0; p < 8; p++) {
        xx += Ni[p] * xn[p];
        yy += Ni[p] * yn[p];
        zz += Ni[p] * zn[p];
      }
      XX[dim * i]     = xx;
      XX[dim * i + 1] = yy;
      XX[dim * i + 2] = zz;
    }

    if (idx == 1) {
      CCmplx zeta, t1, t2;

      xx = XX[dim * i] - 0.8;
      yy = XX[dim * i + 1] + 1.5;

      zeta.real = PetscRealPart(xx);
      zeta.imag = PetscRealPart(yy);

      t1 = CCmplxPow(zeta, -1.0);
      t2 = CCmplxAdd(zeta, t1);

      XX[dim * i]     = CCmplxRe(t2);
      XX[dim * i + 1] = CCmplxIm(t2);
    } else if (idx == 2) {
      CCmplx zeta, t1;

      xx        = XX[dim * i];
      yy        = XX[dim * i + 1];
      zeta.real = PetscRealPart(xx);
      zeta.imag = PetscRealPart(yy);

      t1              = CCmplxSqrt(zeta);
      XX[dim * i]     = CCmplxRe(t1);
      XX[dim * i + 1] = CCmplxIm(t1);
    } else if (idx == 3) {
      CCmplx zeta, t1, t2;

      xx = XX[dim * i] - 0.8;
      yy = XX[dim * i + 1] + 1.5;

      zeta.real       = PetscRealPart(xx);
      zeta.imag       = PetscRealPart(yy);
      t1              = CCmplxPow(zeta, -1.0);
      t2              = CCmplxAdd(zeta, t1);
      XX[dim * i]     = CCmplxRe(t2);
      XX[dim * i + 1] = CCmplxIm(t2);

      xx              = XX[dim * i];
      yy              = XX[dim * i + 1];
      zeta.real       = PetscRealPart(xx);
      zeta.imag       = PetscRealPart(yy);
      t1              = CCmplxExp(zeta);
      XX[dim * i]     = CCmplxRe(t1);
      XX[dim * i + 1] = CCmplxIm(t1);

      xx              = XX[dim * i] + 0.4;
      yy              = XX[dim * i + 1];
      zeta.real       = PetscRealPart(xx);
      zeta.imag       = PetscRealPart(yy);
      t1              = CCmplxPow(zeta, 2.0);
      XX[dim * i]     = CCmplxRe(t1);
      XX[dim * i + 1] = CCmplxIm(t1);
    } else if (idx == 4) {
      PetscScalar Ni[4];
      PetscScalar xi   = XX[dim * i];
      PetscScalar eta  = XX[dim * i + 1];
      PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5};
      PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0};
      PetscInt    p;

      Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
      Ni[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
      Ni[2] = 0.25 * (1.0 - xi) * (1.0 + eta);
      Ni[3] = 0.25 * (1.0 + xi) * (1.0 + eta);

      xx = yy = 0.0;
      for (p = 0; p < 4; p++) {
        xx += Ni[p] * xn[p];
        yy += Ni[p] * yn[p];
      }
      XX[dim * i]     = xx;
      XX[dim * i + 1] = yy;
    }
  }
  PetscCall(VecRestoreArray(Gcoords, &XX));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode DAApplyTrilinearMapping(DM da)
{
  PetscInt      i, j, k;
  PetscInt      sx, nx, sy, ny, sz, nz;
  Vec           Gcoords;
  DMDACoor3d ***XX;
  PetscScalar   xx, yy, zz;
  DM            cda;

  PetscFunctionBeginUser;
  PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
  PetscCall(DMGetCoordinateDM(da, &cda));
  PetscCall(DMGetCoordinates(da, &Gcoords));

  PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
  PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));

  for (i = sx; i < sx + nx; i++) {
    for (j = sy; j < sy + ny; j++) {
      for (k = sz; k < sz + nz; k++) {
        PetscScalar Ni[8];
        PetscScalar xi   = XX[k][j][i].x;
        PetscScalar eta  = XX[k][j][i].y;
        PetscScalar zeta = XX[k][j][i].z;
        PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5, 0.0, 2.1, 0.23, 3.125};
        PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0, -1.45, -0.1, 2.24, 3.79};
        PetscScalar zn[] = {0.0, 0.3, -0.1, 0.123, 0.956, 1.32, 1.12, 0.798};
        PetscInt    p;

        Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
        Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
        Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
        Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);

        Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
        Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
        Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
        Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);

        xx = yy = zz = 0.0;
        for (p = 0; p < 8; p++) {
          xx += Ni[p] * xn[p];
          yy += Ni[p] * yn[p];
          zz += Ni[p] * zn[p];
        }
        XX[k][j][i].x = xx;
        XX[k][j][i].y = yy;
        XX[k][j][i].z = zz;
      }
    }
  }
  PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode DADefineXLinearField2D(DM da, Vec field)
{
  PetscInt      i, j;
  PetscInt      sx, nx, sy, ny;
  Vec           Gcoords;
  DMDACoor2d  **XX;
  PetscScalar **FF;
  DM            cda;

  PetscFunctionBeginUser;
  PetscCall(DMGetCoordinateDM(da, &cda));
  PetscCall(DMGetCoordinates(da, &Gcoords));

  PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
  PetscCall(DMDAVecGetArray(da, field, &FF));

  PetscCall(DMDAGetCorners(da, &sx, &sy, 0, &nx, &ny, 0));

  for (i = sx; i < sx + nx; i++) {
    for (j = sy; j < sy + ny; j++) FF[j][i] = 10.0 + 3.0 * XX[j][i].x + 5.5 * XX[j][i].y + 8.003 * XX[j][i].x * XX[j][i].y;
  }

  PetscCall(DMDAVecRestoreArray(da, field, &FF));
  PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode DADefineXLinearField3D(DM da, Vec field)
{
  PetscInt       i, j, k;
  PetscInt       sx, nx, sy, ny, sz, nz;
  Vec            Gcoords;
  DMDACoor3d  ***XX;
  PetscScalar ***FF;
  DM             cda;

  PetscFunctionBeginUser;
  PetscCall(DMGetCoordinateDM(da, &cda));
  PetscCall(DMGetCoordinates(da, &Gcoords));

  PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
  PetscCall(DMDAVecGetArray(da, field, &FF));

  PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));

  for (k = sz; k < sz + nz; k++) {
    for (j = sy; j < sy + ny; j++) {
      for (i = sx; i < sx + nx; i++) {
        FF[k][j][i] = 10.0 + 4.05 * XX[k][j][i].x + 5.50 * XX[k][j][i].y + 1.33 * XX[k][j][i].z + 2.03 * XX[k][j][i].x * XX[k][j][i].y + 0.03 * XX[k][j][i].x * XX[k][j][i].z + 0.83 * XX[k][j][i].y * XX[k][j][i].z +
                      3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
      }
    }
  }

  PetscCall(DMDAVecRestoreArray(da, field, &FF));
  PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
{
  DM          dac, daf;
  PetscViewer vv;
  Vec         ac, af;
  PetscInt    Mx;
  Mat         II, INTERP;
  Vec         scale;
  PetscBool   output = PETSC_FALSE;

  PetscFunctionBeginUser;
  PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, mx + 1, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, &dac));
  PetscCall(DMSetFromOptions(dac));
  PetscCall(DMSetUp(dac));

  PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
  PetscCall(DMDAGetInfo(daf, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
  Mx--;

  PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
  PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));

  {
    DM  cdaf, cdac;
    Vec coordsc, coordsf;

    PetscCall(DMGetCoordinateDM(dac, &cdac));
    PetscCall(DMGetCoordinateDM(daf, &cdaf));

    PetscCall(DMGetCoordinates(dac, &coordsc));
    PetscCall(DMGetCoordinates(daf, &coordsf));

    PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
    PetscCall(MatInterpolate(II, coordsc, coordsf));
    PetscCall(MatDestroy(&II));
    PetscCall(VecDestroy(&scale));
  }

  PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));

  PetscCall(DMCreateGlobalVector(dac, &ac));
  PetscCall(VecSet(ac, 66.99));

  PetscCall(DMCreateGlobalVector(daf, &af));

  PetscCall(MatMult(INTERP, ac, af));

  {
    Vec       afexact;
    PetscReal nrm;
    PetscInt  N;

    PetscCall(DMCreateGlobalVector(daf, &afexact));
    PetscCall(VecSet(afexact, 66.99));
    PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
    PetscCall(VecNorm(afexact, NORM_2, &nrm));
    PetscCall(VecGetSize(afexact, &N));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "=>%" PetscInt_FMT ", interp err = %1.4e\n", mx, Mx, (double)(nrm / PetscSqrtReal((PetscReal)N))));
    PetscCall(VecDestroy(&afexact));
  }

  PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
  if (output) {
    PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv));
    PetscCall(VecView(ac, vv));
    PetscCall(PetscViewerDestroy(&vv));

    PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv));
    PetscCall(VecView(af, vv));
    PetscCall(PetscViewerDestroy(&vv));
  }

  PetscCall(MatDestroy(&INTERP));
  PetscCall(DMDestroy(&dac));
  PetscCall(DMDestroy(&daf));
  PetscCall(VecDestroy(&ac));
  PetscCall(VecDestroy(&af));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode da_test_RefineCoords2D(PetscInt mx, PetscInt my)
{
  DM          dac, daf;
  PetscViewer vv;
  Vec         ac, af;
  PetscInt    map_id, Mx, My;
  Mat         II, INTERP;
  Vec         scale;
  PetscBool   output = PETSC_FALSE;

  PetscFunctionBeginUser;
  PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, NULL, &dac));
  PetscCall(DMSetFromOptions(dac));
  PetscCall(DMSetUp(dac));

  PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
  PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
  Mx--;
  My--;

  PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
  PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));

  /* apply conformal mappings */
  map_id = 0;
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
  if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));

  {
    DM  cdaf, cdac;
    Vec coordsc, coordsf;

    PetscCall(DMGetCoordinateDM(dac, &cdac));
    PetscCall(DMGetCoordinateDM(daf, &cdaf));

    PetscCall(DMGetCoordinates(dac, &coordsc));
    PetscCall(DMGetCoordinates(daf, &coordsf));

    PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
    PetscCall(MatInterpolate(II, coordsc, coordsf));
    PetscCall(MatDestroy(&II));
    PetscCall(VecDestroy(&scale));
  }

  PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));

  PetscCall(DMCreateGlobalVector(dac, &ac));
  PetscCall(DADefineXLinearField2D(dac, ac));

  PetscCall(DMCreateGlobalVector(daf, &af));
  PetscCall(MatMult(INTERP, ac, af));

  {
    Vec       afexact;
    PetscReal nrm;
    PetscInt  N;

    PetscCall(DMCreateGlobalVector(daf, &afexact));
    PetscCall(VecZeroEntries(afexact));
    PetscCall(DADefineXLinearField2D(daf, afexact));
    PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
    PetscCall(VecNorm(afexact, NORM_2, &nrm));
    PetscCall(VecGetSize(afexact, &N));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT " x %" PetscInt_FMT "]=>[%" PetscInt_FMT " x %" PetscInt_FMT "], interp err = %1.4e\n", mx, my, Mx, My, (double)(nrm / PetscSqrtReal((PetscReal)N))));
    PetscCall(VecDestroy(&afexact));
  }

  PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
  if (output) {
    PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv));
    PetscCall(VecView(ac, vv));
    PetscCall(PetscViewerDestroy(&vv));

    PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv));
    PetscCall(VecView(af, vv));
    PetscCall(PetscViewerDestroy(&vv));
  }

  PetscCall(MatDestroy(&INTERP));
  PetscCall(DMDestroy(&dac));
  PetscCall(DMDestroy(&daf));
  PetscCall(VecDestroy(&ac));
  PetscCall(VecDestroy(&af));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode da_test_RefineCoords3D(PetscInt mx, PetscInt my, PetscInt mz)
{
  DM          dac, daf;
  PetscViewer vv;
  Vec         ac, af;
  PetscInt    map_id, Mx, My, Mz;
  Mat         II, INTERP;
  Vec         scale;
  PetscBool   output = PETSC_FALSE;

  PetscFunctionBeginUser;
  PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, /* 1 dof */
                         1, /* stencil = 1 */ NULL, NULL, NULL, &dac));
  PetscCall(DMSetFromOptions(dac));
  PetscCall(DMSetUp(dac));

  PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
  PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, &Mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
  Mx--;
  My--;
  Mz--;

  PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
  PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));

  /* apply trilinear mappings */
  /*PetscCall(DAApplyTrilinearMapping(dac));*/
  /* apply conformal mappings */
  map_id = 0;
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
  if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));

  {
    DM  cdaf, cdac;
    Vec coordsc, coordsf;

    PetscCall(DMGetCoordinateDM(dac, &cdac));
    PetscCall(DMGetCoordinateDM(daf, &cdaf));

    PetscCall(DMGetCoordinates(dac, &coordsc));
    PetscCall(DMGetCoordinates(daf, &coordsf));

    PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
    PetscCall(MatInterpolate(II, coordsc, coordsf));
    PetscCall(MatDestroy(&II));
    PetscCall(VecDestroy(&scale));
  }

  PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));

  PetscCall(DMCreateGlobalVector(dac, &ac));
  PetscCall(VecZeroEntries(ac));
  PetscCall(DADefineXLinearField3D(dac, ac));

  PetscCall(DMCreateGlobalVector(daf, &af));
  PetscCall(VecZeroEntries(af));

  PetscCall(MatMult(INTERP, ac, af));

  {
    Vec       afexact;
    PetscReal nrm;
    PetscInt  N;

    PetscCall(DMCreateGlobalVector(daf, &afexact));
    PetscCall(VecZeroEntries(afexact));
    PetscCall(DADefineXLinearField3D(daf, afexact));
    PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
    PetscCall(VecNorm(afexact, NORM_2, &nrm));
    PetscCall(VecGetSize(afexact, &N));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "]=>[%" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "], interp err = %1.4e\n", mx, my, mz, Mx, My, Mz, (double)(nrm / PetscSqrtReal((PetscReal)N))));
    PetscCall(VecDestroy(&afexact));
  }

  PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
  if (output) {
    PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv));
    PetscCall(VecView(ac, vv));
    PetscCall(PetscViewerDestroy(&vv));

    PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv));
    PetscCall(VecView(af, vv));
    PetscCall(PetscViewerDestroy(&vv));
  }

  PetscCall(MatDestroy(&INTERP));
  PetscCall(DMDestroy(&dac));
  PetscCall(DMDestroy(&daf));
  PetscCall(VecDestroy(&ac));
  PetscCall(VecDestroy(&af));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char **argv)
{
  PetscInt mx = 2, my = 2, mz = 2, l, nl, dim;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, 0, help));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0));
  nl = 1;
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, 0));
  dim = 2;
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, 0));

  for (l = 0; l < nl; l++) {
    if (dim == 1) PetscCall(da_test_RefineCoords1D(mx));
    else if (dim == 2) PetscCall(da_test_RefineCoords2D(mx, my));
    else if (dim == 3) PetscCall(da_test_RefineCoords3D(mx, my, mz));
    mx = mx * 2;
    my = my * 2;
    mz = mz * 2;
  }
  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

   test:
      suffix: 1d
      args: -mx 10 -nl 6 -dim 1

   test:
      suffix: 2d
      output_file: output/ex36_2d.out
      args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}}

   test:
      suffix: 2dp1
      nsize: 8
      args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
      timeoutfactor: 2

   test:
      suffix: 2dp2
      nsize: 8
      args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
      timeoutfactor: 2

   test:
      suffix: 3d
      args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3

   test:
      suffix: 3dp1
      nsize: 32
      args: -mx 5 -my 5 -mz 5 -nl 3 -dim 3 -cmap 1 -da_refine_x 1 -da_refine_y 3 -da_refine_z 4

TEST*/