Actual source code: ex1.c

  1: static char help[] = "Test point location in DA using DMSwarm\n";

  3: #include <petscdmda.h>
  4: #include <petscdmswarm.h>

  6: PetscErrorCode DMSwarmPrint(DM sw)
  7: {
  8:   PetscReal  *array;
  9:   PetscInt   *pidArray, *cidArray;
 10:   PetscInt    Np, bs;
 11:   PetscMPIInt rank;

 14:   MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank);
 15:   DMSwarmGetLocalSize(sw, &Np);
 16:   DMSwarmGetField(sw, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
 17:   DMSwarmGetField(sw, DMSwarmField_pid, &bs, NULL, (void **)&pidArray);
 18:   DMSwarmGetField(sw, DMSwarmPICField_cellid, &bs, NULL, (void **)&cidArray);
 19:   for (PetscInt p = 0; p < Np; ++p) {
 20:     const PetscReal th = PetscAtan2Real(array[2 * p + 1], array[2 * p]) / PETSC_PI;
 21:     const PetscReal r  = PetscSqrtReal(array[2 * p + 1] * array[2 * p + 1] + array[2 * p] * array[2 * p]);
 22:     PetscPrintf(PETSC_COMM_SELF, "[%d] p %" PetscInt_FMT " (%+1.4f,%+1.4f) r=%1.2f th=%1.3f*pi cellid=%" PetscInt_FMT "\n", rank, pidArray[p], (double)array[2 * p], (double)array[2 * p + 1], (double)r, (double)th, cidArray[p]);
 23:   }
 24:   DMSwarmRestoreField(sw, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
 25:   DMSwarmRestoreField(sw, DMSwarmField_pid, &bs, NULL, (void **)&pidArray);
 26:   DMSwarmRestoreField(sw, DMSwarmPICField_cellid, &bs, NULL, (void **)&pidArray);
 27:   return 0;
 28: }

 30: int main(int argc, char **argv)
 31: {
 32:   DM            dm, sw;
 33:   PetscDataType dtype;
 34:   PetscReal    *coords, r, dr;
 35:   PetscInt      Nx = 4, Ny = 4, Np = 8, bs;
 36:   PetscMPIInt   rank, size;

 38:   PetscInitialize(&argc, &argv, NULL, help);
 39:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 40:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 42:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, Nx, Ny, PETSC_DECIDE, PETSC_DECIDE, 1, 2, NULL, NULL, &dm);
 43:   DMSetFromOptions(dm);
 44:   DMSetUp(dm);
 45:   DMDASetUniformCoordinates(dm, -1., 1., -1., 1., -1., 1.);
 46:   DMViewFromOptions(dm, NULL, "-da_view");

 48:   DMCreate(PETSC_COMM_WORLD, &sw);
 49:   PetscObjectSetName((PetscObject)sw, "Particle Grid");
 50:   DMSetType(sw, DMSWARM);
 51:   DMSetDimension(sw, 2);
 52:   DMSwarmSetType(sw, DMSWARM_PIC);
 53:   DMSetFromOptions(sw);
 54:   DMSwarmSetCellDM(sw, dm);
 55:   DMSwarmInitializeFieldRegister(sw);
 56:   DMSwarmRegisterPetscDatatypeField(sw, "u", 1, PETSC_SCALAR);
 57:   DMSwarmFinalizeFieldRegister(sw);
 58:   DMSwarmSetLocalSizes(sw, Np, 2);
 59:   DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords);
 60:   dr = 1.0 / (size + 1);
 61:   r  = (rank + 1) * dr;
 62:   for (PetscInt p = 0; p < Np; ++p) {
 63:     const PetscReal th = (p + 0.5) * 2. * PETSC_PI / Np;

 65:     coords[p * 2 + 0] = r * PetscCosReal(th);
 66:     coords[p * 2 + 1] = r * PetscSinReal(th);
 67:   }
 68:   DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords);
 69:   DMViewFromOptions(sw, NULL, "-swarm_view");
 70:   DMSwarmPrint(sw);

 72:   PetscPrintf(PETSC_COMM_WORLD, "\n... calling DMSwarmMigrate ...\n\n");
 73:   DMSwarmMigrate(sw, PETSC_TRUE);
 74:   DMViewFromOptions(sw, NULL, "-swarm_view");
 75:   DMSwarmPrint(sw);

 77:   DMDestroy(&sw);
 78:   DMDestroy(&dm);
 79:   PetscFinalize();
 80:   return 0;
 81: }

 83: /*TEST

 85:   # Swarm does not handle complex or quad
 86:   build:
 87:     requires: !complex double

 89:   test:
 90:     suffix: 0

 92: TEST*/