Actual source code: ex1.c
1: static char help[] = "Test point location in DA using DMSwarm\n";
3: #include <petscdmda.h>
4: #include <petscdmswarm.h>
6: PetscErrorCode DMSwarmPrint(DM sw)
7: {
8: PetscReal *array;
9: PetscInt *pidArray, *cidArray;
10: PetscInt Np, bs;
11: PetscMPIInt rank;
14: MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank);
15: DMSwarmGetLocalSize(sw, &Np);
16: DMSwarmGetField(sw, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
17: DMSwarmGetField(sw, DMSwarmField_pid, &bs, NULL, (void **)&pidArray);
18: DMSwarmGetField(sw, DMSwarmPICField_cellid, &bs, NULL, (void **)&cidArray);
19: for (PetscInt p = 0; p < Np; ++p) {
20: const PetscReal th = PetscAtan2Real(array[2 * p + 1], array[2 * p]) / PETSC_PI;
21: const PetscReal r = PetscSqrtReal(array[2 * p + 1] * array[2 * p + 1] + array[2 * p] * array[2 * p]);
22: PetscPrintf(PETSC_COMM_SELF, "[%d] p %" PetscInt_FMT " (%+1.4f,%+1.4f) r=%1.2f th=%1.3f*pi cellid=%" PetscInt_FMT "\n", rank, pidArray[p], (double)array[2 * p], (double)array[2 * p + 1], (double)r, (double)th, cidArray[p]);
23: }
24: DMSwarmRestoreField(sw, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
25: DMSwarmRestoreField(sw, DMSwarmField_pid, &bs, NULL, (void **)&pidArray);
26: DMSwarmRestoreField(sw, DMSwarmPICField_cellid, &bs, NULL, (void **)&pidArray);
27: return 0;
28: }
30: int main(int argc, char **argv)
31: {
32: DM dm, sw;
33: PetscDataType dtype;
34: PetscReal *coords, r, dr;
35: PetscInt Nx = 4, Ny = 4, Np = 8, bs;
36: PetscMPIInt rank, size;
38: PetscInitialize(&argc, &argv, NULL, help);
39: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
40: MPI_Comm_size(PETSC_COMM_WORLD, &size);
42: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, Nx, Ny, PETSC_DECIDE, PETSC_DECIDE, 1, 2, NULL, NULL, &dm);
43: DMSetFromOptions(dm);
44: DMSetUp(dm);
45: DMDASetUniformCoordinates(dm, -1., 1., -1., 1., -1., 1.);
46: DMViewFromOptions(dm, NULL, "-da_view");
48: DMCreate(PETSC_COMM_WORLD, &sw);
49: PetscObjectSetName((PetscObject)sw, "Particle Grid");
50: DMSetType(sw, DMSWARM);
51: DMSetDimension(sw, 2);
52: DMSwarmSetType(sw, DMSWARM_PIC);
53: DMSetFromOptions(sw);
54: DMSwarmSetCellDM(sw, dm);
55: DMSwarmInitializeFieldRegister(sw);
56: DMSwarmRegisterPetscDatatypeField(sw, "u", 1, PETSC_SCALAR);
57: DMSwarmFinalizeFieldRegister(sw);
58: DMSwarmSetLocalSizes(sw, Np, 2);
59: DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords);
60: dr = 1.0 / (size + 1);
61: r = (rank + 1) * dr;
62: for (PetscInt p = 0; p < Np; ++p) {
63: const PetscReal th = (p + 0.5) * 2. * PETSC_PI / Np;
65: coords[p * 2 + 0] = r * PetscCosReal(th);
66: coords[p * 2 + 1] = r * PetscSinReal(th);
67: }
68: DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords);
69: DMViewFromOptions(sw, NULL, "-swarm_view");
70: DMSwarmPrint(sw);
72: PetscPrintf(PETSC_COMM_WORLD, "\n... calling DMSwarmMigrate ...\n\n");
73: DMSwarmMigrate(sw, PETSC_TRUE);
74: DMViewFromOptions(sw, NULL, "-swarm_view");
75: DMSwarmPrint(sw);
77: DMDestroy(&sw);
78: DMDestroy(&dm);
79: PetscFinalize();
80: return 0;
81: }
83: /*TEST
85: # Swarm does not handle complex or quad
86: build:
87: requires: !complex double
89: test:
90: suffix: 0
92: TEST*/