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*/