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";
2
3 #include <petscdmplex.h>
4 #include <petscdmswarm.h>
5 #include <petscts.h>
6
7 typedef struct {
8 PetscInt particlesPerCell; /* The number of partices per cell */
9 } AppCtx;
10
ProcessOptions(MPI_Comm comm,AppCtx * options)11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
12 {
13 PetscFunctionBeginUser;
14 options->particlesPerCell = 1;
15
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 }
21
CreateMesh(MPI_Comm comm,DM * dm,AppCtx * user)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 }
31
CreateParticles(DM dm,DM * sw,AppCtx * user)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;
38
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;
57
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 }
66
main(int argc,char ** argv)67 int main(int argc, char **argv)
68 {
69 DM dm, sw, cellsw; /* Mesh and particle managers */
70 MPI_Comm comm;
71 AppCtx user;
72
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 }
91
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*/
99