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:   PetscInt *cellid;
 35:   PetscInt  dim, cStart, cEnd, c, Np = user->particlesPerCell, p;

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

 54:       cellid[n] = c;
 55:     }
 56:   }
 57:   PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
 58:   PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
 59:   PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
 60:   PetscFunctionReturn(PETSC_SUCCESS);
 61: }

 63: int main(int argc, char **argv)
 64: {
 65:   DM       dm, sw, cellsw; /* Mesh and particle managers */
 66:   MPI_Comm comm;
 67:   AppCtx   user;

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

 88: /*TEST
 89:    build:
 90:      requires: triangle !single !complex
 91:    test:
 92:      suffix: 1
 93:      args: -particles_per_cell 2 -dm_plex_box_faces 2,1 -dm_view -sw_view -subswarm_view
 94: TEST*/