Actual source code: ex3.c

  1: static char help[] = "Example usage of extracting single cells with their associated fields from a swarm and putting it in a new swarm object\n";

  3: #include <petscdmplex.h>
  4: #include <petscdmswarm.h>
  5: #include <petscts.h>

  7: typedef struct {
  8:   PetscInt particlesPerCell; /* The number of partices per cell */
  9: } AppCtx;

 11: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 12: {
 13:   PetscFunctionBeginUser;
 14:   options->particlesPerCell = 1;

 16:   PetscOptionsBegin(comm, "", "CellSwarm Options", "DMSWARM");
 17:   PetscCall(PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex3.c", options->particlesPerCell, &options->particlesPerCell, NULL));
 18:   PetscOptionsEnd();
 19:   PetscFunctionReturn(PETSC_SUCCESS);
 20: }

 22: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
 23: {
 24:   PetscFunctionBeginUser;
 25:   PetscCall(DMCreate(comm, dm));
 26:   PetscCall(DMSetType(*dm, DMPLEX));
 27:   PetscCall(DMSetFromOptions(*dm));
 28:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 29:   PetscFunctionReturn(PETSC_SUCCESS);
 30: }

 32: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
 33: {
 34:   DMSwarmCellDM celldm;
 35:   PetscInt     *swarm_cellid;
 36:   PetscInt      dim, cStart, cEnd, c, Np = user->particlesPerCell, p;
 37:   const char   *cellid;

 39:   PetscFunctionBeginUser;
 40:   PetscCall(DMGetDimension(dm, &dim));
 41:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
 42:   PetscCall(DMSetType(*sw, DMSWARM));
 43:   PetscCall(DMSetDimension(*sw, dim));
 44:   PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
 45:   PetscCall(DMSwarmSetCellDM(*sw, dm));
 46:   PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL));
 47:   PetscCall(DMSwarmFinalizeFieldRegister(*sw));
 48:   PetscCall(DMSwarmGetCellDMActive(*sw, &celldm));
 49:   PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
 50:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 51:   PetscCall(DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0));
 52:   PetscCall(DMSetFromOptions(*sw));
 53:   PetscCall(DMSwarmGetField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
 54:   for (c = cStart; c < cEnd; ++c) {
 55:     for (p = 0; p < Np; ++p) {
 56:       const PetscInt n = c * Np + p;

 58:       swarm_cellid[n] = c;
 59:     }
 60:   }
 61:   PetscCall(DMSwarmRestoreField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
 62:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
 63:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
 64:   PetscFunctionReturn(PETSC_SUCCESS);
 65: }

 67: int main(int argc, char **argv)
 68: {
 69:   DM       dm, sw, cellsw; /* Mesh and particle managers */
 70:   MPI_Comm comm;
 71:   AppCtx   user;

 73:   PetscFunctionBeginUser;
 74:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 75:   comm = PETSC_COMM_WORLD;
 76:   PetscCall(ProcessOptions(comm, &user));
 77:   PetscCall(CreateMesh(comm, &dm, &user));
 78:   PetscCall(CreateParticles(dm, &sw, &user));
 79:   PetscCall(DMSetApplicationContext(sw, &user));
 80:   PetscCall(DMCreate(comm, &cellsw));
 81:   PetscCall(PetscObjectSetName((PetscObject)cellsw, "SubParticles"));
 82:   PetscCall(DMSwarmGetCellSwarm(sw, 1, cellsw));
 83:   PetscCall(DMViewFromOptions(cellsw, NULL, "-subswarm_view"));
 84:   PetscCall(DMSwarmRestoreCellSwarm(sw, 1, cellsw));
 85:   PetscCall(DMDestroy(&sw));
 86:   PetscCall(DMDestroy(&dm));
 87:   PetscCall(DMDestroy(&cellsw));
 88:   PetscCall(PetscFinalize());
 89:   return 0;
 90: }

 92: /*TEST
 93:    build:
 94:      requires: triangle !single !complex
 95:    test:
 96:      suffix: 1
 97:      args: -particles_per_cell 2 -dm_plex_box_faces 2,1 -dm_view -sw_view -subswarm_view
 98: TEST*/