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