static char help[] = "Tests DMSwarm\n\n";

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

/*
 Checks for variable blocksize
*/
PetscErrorCode ex2_1(void)
{
  DM          dms;
  Vec         x;
  PetscMPIInt rank;
  PetscInt    p, bs, nlocal;

  PetscFunctionBegin;
  PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

  PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
  PetscCall(DMSetType(dms, DMSWARM));
  PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));
  PetscCall(DMSwarmInitializeFieldRegister(dms));
  PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
  PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "strain", 3, PETSC_REAL));
  PetscCall(DMSwarmFinalizeFieldRegister(dms));
  PetscCall(DMSwarmSetLocalSizes(dms, 5 + rank, 4));
  PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(DMSwarmGetLocalSize(dms, &nlocal));

  {
    PetscReal *array;
    PetscCall(DMSwarmGetField(dms, "viscosity", &bs, NULL, (void **)&array));
    for (p = 0; p < nlocal; p++) array[p] = 11.1 + p * 0.1 + rank * 100.0;
    PetscCall(DMSwarmRestoreField(dms, "viscosity", &bs, NULL, (void **)&array));
  }

  {
    PetscReal *array;
    PetscCall(DMSwarmGetField(dms, "strain", &bs, NULL, (void **)&array));
    for (p = 0; p < nlocal; p++) {
      array[bs * p + 0] = 2.0e-2 + p * 0.001 + rank * 1.0;
      array[bs * p + 1] = 2.0e-2 + p * 0.002 + rank * 1.0;
      array[bs * p + 2] = 2.0e-2 + p * 0.003 + rank * 1.0;
    }
    PetscCall(DMSwarmRestoreField(dms, "strain", &bs, NULL, (void **)&array));
  }

  PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
  PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));

  PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "strain", &x));
  PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "strain", &x));

  PetscCall(DMSwarmVectorDefineField(dms, "strain"));
  PetscCall(DMCreateGlobalVector(dms, &x));
  PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
  PetscCall(VecDestroy(&x));
  PetscCall(DMDestroy(&dms));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

/*TEST

   test:
      requires: !complex double
      nsize: 3
      filter: grep -v atomic
      filter_output: grep -v atomic

TEST*/