Actual source code: swarm_ex1.c

  1: static char help[] = "Tests DMSwarm\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmda.h>
  5: #include <petscdmswarm.h>

  7: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM, PetscErrorCode (*)(DM, void *, PetscInt *, PetscInt **), size_t, void *, PetscInt *);
  8: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM, PetscInt *);

 10: PetscErrorCode ex1_1(void)
 11: {
 12:   DM          dms;
 13:   Vec         x;
 14:   PetscMPIInt rank, size;
 15:   PetscInt    p;

 17:   PetscFunctionBegin;
 18:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 19:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 20:   PetscCheck(!(size > 1) || !(size != 4), PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must be run with 4 MPI ranks");

 22:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
 23:   PetscCall(DMSetType(dms, DMSWARM));
 24:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));

 26:   PetscCall(DMSwarmInitializeFieldRegister(dms));
 27:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
 28:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "strain", 1, PETSC_REAL));
 29:   PetscCall(DMSwarmFinalizeFieldRegister(dms));
 30:   PetscCall(DMSwarmSetLocalSizes(dms, 5 + rank, 4));
 31:   PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));

 33:   {
 34:     PetscReal *array;
 35:     PetscCall(DMSwarmGetField(dms, "viscosity", NULL, NULL, (void **)&array));
 36:     for (p = 0; p < 5 + rank; p++) array[p] = 11.1 + p * 0.1 + rank * 100.0;
 37:     PetscCall(DMSwarmRestoreField(dms, "viscosity", NULL, NULL, (void **)&array));
 38:   }

 40:   {
 41:     PetscReal *array;
 42:     PetscCall(DMSwarmGetField(dms, "strain", NULL, NULL, (void **)&array));
 43:     for (p = 0; p < 5 + rank; p++) array[p] = 2.0e-2 + p * 0.001 + rank * 1.0;
 44:     PetscCall(DMSwarmRestoreField(dms, "strain", NULL, NULL, (void **)&array));
 45:   }

 47:   PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
 48:   PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));

 50:   PetscCall(DMSwarmVectorDefineField(dms, "strain"));
 51:   PetscCall(DMCreateGlobalVector(dms, &x));
 52:   PetscCall(VecDestroy(&x));

 54:   {
 55:     PetscInt *rankval;
 56:     PetscInt  npoints[2], npoints_orig[2];

 58:     PetscCall(DMSwarmGetLocalSize(dms, &npoints_orig[0]));
 59:     PetscCall(DMSwarmGetSize(dms, &npoints_orig[1]));
 60:     PetscCall(DMSwarmGetField(dms, "DMSwarm_rank", NULL, NULL, (void **)&rankval));
 61:     if ((rank == 0) && (size > 1)) {
 62:       rankval[0] = 1;
 63:       rankval[3] = 1;
 64:     }
 65:     if (rank == 3) rankval[2] = 1;
 66:     PetscCall(DMSwarmRestoreField(dms, "DMSwarm_rank", NULL, NULL, (void **)&rankval));
 67:     PetscCall(DMSwarmMigrate(dms, PETSC_TRUE));
 68:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
 69:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
 70:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
 71:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ") after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints_orig[0], npoints_orig[1], npoints[0], npoints[1]));
 72:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
 73:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
 74:   }
 75:   {
 76:     PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
 77:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
 78:     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));
 79:   }
 80:   {
 81:     PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "strain", &x));
 82:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
 83:     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "strain", &x));
 84:   }

 86:   PetscCall(DMDestroy(&dms));
 87:   PetscFunctionReturn(PETSC_SUCCESS);
 88: }

 90: PetscErrorCode ex1_2(void)
 91: {
 92:   DM          dms;
 93:   Vec         x;
 94:   PetscMPIInt rank, size;
 95:   PetscInt    p;

 97:   PetscFunctionBegin;
 98:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 99:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
100:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
101:   PetscCall(DMSetType(dms, DMSWARM));
102:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));
103:   PetscCall(DMSwarmInitializeFieldRegister(dms));

105:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
106:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "strain", 1, PETSC_REAL));
107:   PetscCall(DMSwarmFinalizeFieldRegister(dms));
108:   PetscCall(DMSwarmSetLocalSizes(dms, 5 + rank, 4));
109:   PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));
110:   {
111:     PetscReal *array;
112:     PetscCall(DMSwarmGetField(dms, "viscosity", NULL, NULL, (void **)&array));
113:     for (p = 0; p < 5 + rank; p++) array[p] = 11.1 + p * 0.1 + rank * 100.0;
114:     PetscCall(DMSwarmRestoreField(dms, "viscosity", NULL, NULL, (void **)&array));
115:   }
116:   {
117:     PetscReal *array;
118:     PetscCall(DMSwarmGetField(dms, "strain", NULL, NULL, (void **)&array));
119:     for (p = 0; p < 5 + rank; p++) array[p] = 2.0e-2 + p * 0.001 + rank * 1.0;
120:     PetscCall(DMSwarmRestoreField(dms, "strain", NULL, NULL, (void **)&array));
121:   }
122:   {
123:     PetscInt *rankval;
124:     PetscInt  npoints[2], npoints_orig[2];

126:     PetscCall(DMSwarmGetLocalSize(dms, &npoints_orig[0]));
127:     PetscCall(DMSwarmGetSize(dms, &npoints_orig[1]));
128:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
129:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints_orig[0], npoints_orig[1]));
130:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
131:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));

133:     PetscCall(DMSwarmGetField(dms, "DMSwarm_rank", NULL, NULL, (void **)&rankval));

135:     if (rank == 1) rankval[0] = -1;
136:     if (rank == 2) rankval[1] = -1;
137:     if (rank == 3) {
138:       rankval[3] = -1;
139:       rankval[4] = -1;
140:     }
141:     PetscCall(DMSwarmRestoreField(dms, "DMSwarm_rank", NULL, NULL, (void **)&rankval));
142:     PetscCall(DMSwarmCollectViewCreate(dms));
143:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
144:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
145:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
146:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints[0], npoints[1]));
147:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
148:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));

150:     PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
151:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
152:     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));

154:     PetscCall(DMSwarmCollectViewDestroy(dms));
155:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
156:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
157:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
158:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] after_v(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints[0], npoints[1]));
159:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
160:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));

162:     PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
163:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
164:     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));
165:   }
166:   PetscCall(DMDestroy(&dms));
167:   PetscFunctionReturn(PETSC_SUCCESS);
168: }

170: /*
171:  splot "c-rank0.gp","c-rank1.gp","c-rank2.gp","c-rank3.gp"
172: */
173: PetscErrorCode ex1_3(void)
174: {
175:   DM          dms;
176:   PetscMPIInt rank, size;
177:   PetscInt    is, js, ni, nj, overlap;
178:   DM          dmcell;

180:   PetscFunctionBegin;
181:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
182:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
183:   overlap = 2;
184:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 13, 13, PETSC_DECIDE, PETSC_DECIDE, 1, overlap, NULL, NULL, &dmcell));
185:   PetscCall(DMSetFromOptions(dmcell));
186:   PetscCall(DMSetUp(dmcell));
187:   PetscCall(DMDASetUniformCoordinates(dmcell, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
188:   PetscCall(DMDAGetCorners(dmcell, &is, &js, NULL, &ni, &nj, NULL));
189:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
190:   PetscCall(DMSetType(dms, DMSWARM));
191:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));
192:   PetscCall(DMSwarmSetCellDM(dms, dmcell));

194:   /* load in data types */
195:   PetscCall(DMSwarmInitializeFieldRegister(dms));
196:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
197:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "coorx", 1, PETSC_REAL));
198:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "coory", 1, PETSC_REAL));
199:   PetscCall(DMSwarmFinalizeFieldRegister(dms));
200:   PetscCall(DMSwarmSetLocalSizes(dms, ni * nj * 4, 4));
201:   PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));

203:   /* set values within the swarm */
204:   {
205:     PetscReal   *array_x, *array_y;
206:     PetscInt     npoints, i, j, cnt;
207:     DMDACoor2d **LA_coor;
208:     Vec          coor;
209:     DM           dmcellcdm;

211:     PetscCall(DMGetCoordinateDM(dmcell, &dmcellcdm));
212:     PetscCall(DMGetCoordinates(dmcell, &coor));
213:     PetscCall(DMDAVecGetArray(dmcellcdm, coor, &LA_coor));
214:     PetscCall(DMSwarmGetLocalSize(dms, &npoints));
215:     PetscCall(DMSwarmGetField(dms, "coorx", NULL, NULL, (void **)&array_x));
216:     PetscCall(DMSwarmGetField(dms, "coory", NULL, NULL, (void **)&array_y));
217:     cnt = 0;
218:     for (j = js; j < js + nj; j++) {
219:       for (i = is; i < is + ni; i++) {
220:         PetscReal xp, yp;
221:         xp                   = PetscRealPart(LA_coor[j][i].x);
222:         yp                   = PetscRealPart(LA_coor[j][i].y);
223:         array_x[4 * cnt + 0] = xp - 0.05;
224:         if (array_x[4 * cnt + 0] < -1.0) array_x[4 * cnt + 0] = -1.0 + 1.0e-12;
225:         array_x[4 * cnt + 1] = xp + 0.05;
226:         if (array_x[4 * cnt + 1] > 1.0) array_x[4 * cnt + 1] = 1.0 - 1.0e-12;
227:         array_x[4 * cnt + 2] = xp - 0.05;
228:         if (array_x[4 * cnt + 2] < -1.0) array_x[4 * cnt + 2] = -1.0 + 1.0e-12;
229:         array_x[4 * cnt + 3] = xp + 0.05;
230:         if (array_x[4 * cnt + 3] > 1.0) array_x[4 * cnt + 3] = 1.0 - 1.0e-12;

232:         array_y[4 * cnt + 0] = yp - 0.05;
233:         if (array_y[4 * cnt + 0] < -1.0) array_y[4 * cnt + 0] = -1.0 + 1.0e-12;
234:         array_y[4 * cnt + 1] = yp - 0.05;
235:         if (array_y[4 * cnt + 1] < -1.0) array_y[4 * cnt + 1] = -1.0 + 1.0e-12;
236:         array_y[4 * cnt + 2] = yp + 0.05;
237:         if (array_y[4 * cnt + 2] > 1.0) array_y[4 * cnt + 2] = 1.0 - 1.0e-12;
238:         array_y[4 * cnt + 3] = yp + 0.05;
239:         if (array_y[4 * cnt + 3] > 1.0) array_y[4 * cnt + 3] = 1.0 - 1.0e-12;
240:         cnt++;
241:       }
242:     }
243:     PetscCall(DMSwarmRestoreField(dms, "coory", NULL, NULL, (void **)&array_y));
244:     PetscCall(DMSwarmRestoreField(dms, "coorx", NULL, NULL, (void **)&array_x));
245:     PetscCall(DMDAVecRestoreArray(dmcellcdm, coor, &LA_coor));
246:   }
247:   {
248:     PetscInt npoints[2], npoints_orig[2], ng;

250:     PetscCall(DMSwarmGetLocalSize(dms, &npoints_orig[0]));
251:     PetscCall(DMSwarmGetSize(dms, &npoints_orig[1]));
252:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
253:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints_orig[0], npoints_orig[1]));
254:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
255:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
256:     PetscCall(DMSwarmCollect_DMDABoundingBox(dms, &ng));

258:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
259:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
260:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
261:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints[0], npoints[1]));
262:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
263:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
264:   }
265:   {
266:     PetscReal *array_x, *array_y;
267:     PetscInt   npoints, p;
268:     FILE      *fp = NULL;
269:     char       name[PETSC_MAX_PATH_LEN];

271:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "c-rank%d.gp", rank));
272:     fp = fopen(name, "w");
273:     PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file %s", name);
274:     PetscCall(DMSwarmGetLocalSize(dms, &npoints));
275:     PetscCall(DMSwarmGetField(dms, "coorx", NULL, NULL, (void **)&array_x));
276:     PetscCall(DMSwarmGetField(dms, "coory", NULL, NULL, (void **)&array_y));
277:     for (p = 0; p < npoints; p++) fprintf(fp, "%+1.4e %+1.4e %1.4e\n", array_x[p], array_y[p], (double)rank);
278:     PetscCall(DMSwarmRestoreField(dms, "coory", NULL, NULL, (void **)&array_y));
279:     PetscCall(DMSwarmRestoreField(dms, "coorx", NULL, NULL, (void **)&array_x));
280:     fclose(fp);
281:   }
282:   PetscCall(DMDestroy(&dmcell));
283:   PetscCall(DMDestroy(&dms));
284:   PetscFunctionReturn(PETSC_SUCCESS);
285: }

287: typedef struct {
288:   PetscReal cx[2];
289:   PetscReal radius;
290: } CollectZoneCtx;

292: PetscErrorCode collect_zone(DM dm, void *ctx, PetscInt *nfound, PetscInt **foundlist)
293: {
294:   CollectZoneCtx *zone = (CollectZoneCtx *)ctx;
295:   PetscInt        p, npoints;
296:   PetscReal      *array_x, *array_y, r2;
297:   PetscInt        p2collect, *plist;
298:   PetscMPIInt     rank;

300:   PetscFunctionBegin;
301:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
302:   PetscCall(DMSwarmGetLocalSize(dm, &npoints));
303:   PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
304:   PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));

306:   r2        = zone->radius * zone->radius;
307:   p2collect = 0;
308:   for (p = 0; p < npoints; p++) {
309:     PetscReal sep2;

311:     sep2 = (array_x[p] - zone->cx[0]) * (array_x[p] - zone->cx[0]);
312:     sep2 += (array_y[p] - zone->cx[1]) * (array_y[p] - zone->cx[1]);
313:     if (sep2 < r2) p2collect++;
314:   }

316:   PetscCall(PetscMalloc1(p2collect + 1, &plist));
317:   p2collect = 0;
318:   for (p = 0; p < npoints; p++) {
319:     PetscReal sep2;

321:     sep2 = (array_x[p] - zone->cx[0]) * (array_x[p] - zone->cx[0]);
322:     sep2 += (array_y[p] - zone->cx[1]) * (array_y[p] - zone->cx[1]);
323:     if (sep2 < r2) {
324:       plist[p2collect] = p;
325:       p2collect++;
326:     }
327:   }
328:   PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
329:   PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));

331:   *nfound    = p2collect;
332:   *foundlist = plist;
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

336: PetscErrorCode ex1_4(void)
337: {
338:   DM              dms;
339:   PetscMPIInt     rank, size;
340:   PetscInt        is, js, ni, nj, overlap, nn;
341:   DM              dmcell;
342:   CollectZoneCtx *zone;
343:   PetscReal       dx;

345:   PetscFunctionBegin;
346:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
347:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
348:   nn      = 101;
349:   dx      = 2.0 / (PetscReal)(nn - 1);
350:   overlap = 0;
351:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, nn, nn, PETSC_DECIDE, PETSC_DECIDE, 1, overlap, NULL, NULL, &dmcell));
352:   PetscCall(DMSetFromOptions(dmcell));
353:   PetscCall(DMSetUp(dmcell));
354:   PetscCall(DMDASetUniformCoordinates(dmcell, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
355:   PetscCall(DMDAGetCorners(dmcell, &is, &js, NULL, &ni, &nj, NULL));
356:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
357:   PetscCall(DMSetType(dms, DMSWARM));
358:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));

360:   /* load in data types */
361:   PetscCall(DMSwarmInitializeFieldRegister(dms));
362:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
363:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "coorx", 1, PETSC_REAL));
364:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "coory", 1, PETSC_REAL));
365:   PetscCall(DMSwarmFinalizeFieldRegister(dms));
366:   PetscCall(DMSwarmSetLocalSizes(dms, ni * nj * 4, 4));
367:   PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));

369:   /* set values within the swarm */
370:   {
371:     PetscReal   *array_x, *array_y;
372:     PetscInt     npoints, i, j, cnt;
373:     DMDACoor2d **LA_coor;
374:     Vec          coor;
375:     DM           dmcellcdm;

377:     PetscCall(DMGetCoordinateDM(dmcell, &dmcellcdm));
378:     PetscCall(DMGetCoordinates(dmcell, &coor));
379:     PetscCall(DMDAVecGetArray(dmcellcdm, coor, &LA_coor));
380:     PetscCall(DMSwarmGetLocalSize(dms, &npoints));
381:     PetscCall(DMSwarmGetField(dms, "coorx", NULL, NULL, (void **)&array_x));
382:     PetscCall(DMSwarmGetField(dms, "coory", NULL, NULL, (void **)&array_y));
383:     cnt = 0;
384:     for (j = js; j < js + nj; j++) {
385:       for (i = is; i < is + ni; i++) {
386:         PetscReal xp, yp;

388:         xp                   = PetscRealPart(LA_coor[j][i].x);
389:         yp                   = PetscRealPart(LA_coor[j][i].y);
390:         array_x[4 * cnt + 0] = xp - dx * 0.1; /*if (array_x[4*cnt+0] < -1.0) array_x[4*cnt+0] = -1.0+1.0e-12;*/
391:         array_x[4 * cnt + 1] = xp + dx * 0.1; /*if (array_x[4*cnt+1] > 1.0)  { array_x[4*cnt+1] =  1.0-1.0e-12; }*/
392:         array_x[4 * cnt + 2] = xp - dx * 0.1; /*if (array_x[4*cnt+2] < -1.0) array_x[4*cnt+2] = -1.0+1.0e-12;*/
393:         array_x[4 * cnt + 3] = xp + dx * 0.1; /*if (array_x[4*cnt+3] > 1.0)  { array_x[4*cnt+3] =  1.0-1.0e-12; }*/
394:         array_y[4 * cnt + 0] = yp - dx * 0.1; /*if (array_y[4*cnt+0] < -1.0) array_y[4*cnt+0] = -1.0+1.0e-12;*/
395:         array_y[4 * cnt + 1] = yp - dx * 0.1; /*if (array_y[4*cnt+1] < -1.0) array_y[4*cnt+1] = -1.0+1.0e-12;*/
396:         array_y[4 * cnt + 2] = yp + dx * 0.1; /*if (array_y[4*cnt+2] > 1.0)  { array_y[4*cnt+2] =  1.0-1.0e-12; }*/
397:         array_y[4 * cnt + 3] = yp + dx * 0.1; /*if (array_y[4*cnt+3] > 1.0)  { array_y[4*cnt+3] =  1.0-1.0e-12; }*/
398:         cnt++;
399:       }
400:     }
401:     PetscCall(DMSwarmRestoreField(dms, "coory", NULL, NULL, (void **)&array_y));
402:     PetscCall(DMSwarmRestoreField(dms, "coorx", NULL, NULL, (void **)&array_x));
403:     PetscCall(DMDAVecRestoreArray(dmcellcdm, coor, &LA_coor));
404:   }
405:   PetscCall(PetscMalloc1(1, &zone));
406:   if (size == 4) {
407:     if (rank == 0) {
408:       zone->cx[0]  = 0.5;
409:       zone->cx[1]  = 0.5;
410:       zone->radius = 0.3;
411:     }
412:     if (rank == 1) {
413:       zone->cx[0]  = -0.5;
414:       zone->cx[1]  = 0.5;
415:       zone->radius = 0.25;
416:     }
417:     if (rank == 2) {
418:       zone->cx[0]  = 0.5;
419:       zone->cx[1]  = -0.5;
420:       zone->radius = 0.2;
421:     }
422:     if (rank == 3) {
423:       zone->cx[0]  = -0.5;
424:       zone->cx[1]  = -0.5;
425:       zone->radius = 0.1;
426:     }
427:   } else {
428:     if (rank == 0) {
429:       zone->cx[0]  = 0.5;
430:       zone->cx[1]  = 0.5;
431:       zone->radius = 0.8;
432:     } else {
433:       zone->cx[0]  = 10.0;
434:       zone->cx[1]  = 10.0;
435:       zone->radius = 0.0;
436:     }
437:   }
438:   {
439:     PetscInt npoints[2], npoints_orig[2], ng;

441:     PetscCall(DMSwarmGetLocalSize(dms, &npoints_orig[0]));
442:     PetscCall(DMSwarmGetSize(dms, &npoints_orig[1]));
443:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
444:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints_orig[0], npoints_orig[1]));
445:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
446:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
447:     PetscCall(DMSwarmCollect_General(dms, collect_zone, sizeof(CollectZoneCtx), zone, &ng));
448:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
449:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
450:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
451:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints[0], npoints[1]));
452:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
453:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
454:   }
455:   {
456:     PetscReal *array_x, *array_y;
457:     PetscInt   npoints, p;
458:     FILE      *fp = NULL;
459:     char       name[PETSC_MAX_PATH_LEN];

461:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "c-rank%d.gp", rank));
462:     fp = fopen(name, "w");
463:     PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file %s", name);
464:     PetscCall(DMSwarmGetLocalSize(dms, &npoints));
465:     PetscCall(DMSwarmGetField(dms, "coorx", NULL, NULL, (void **)&array_x));
466:     PetscCall(DMSwarmGetField(dms, "coory", NULL, NULL, (void **)&array_y));
467:     for (p = 0; p < npoints; p++) fprintf(fp, "%+1.4e %+1.4e %1.4e\n", array_x[p], array_y[p], (double)rank);
468:     PetscCall(DMSwarmRestoreField(dms, "coory", NULL, NULL, (void **)&array_y));
469:     PetscCall(DMSwarmRestoreField(dms, "coorx", NULL, NULL, (void **)&array_x));
470:     fclose(fp);
471:   }
472:   PetscCall(DMDestroy(&dmcell));
473:   PetscCall(DMDestroy(&dms));
474:   PetscCall(PetscFree(zone));
475:   PetscFunctionReturn(PETSC_SUCCESS);
476: }

478: int main(int argc, char **argv)
479: {
480:   PetscInt test_mode = 4;

482:   PetscFunctionBeginUser;
483:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
484:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_mode", &test_mode, NULL));
485:   if (test_mode == 1) {
486:     PetscCall(ex1_1());
487:   } else if (test_mode == 2) {
488:     PetscCall(ex1_2());
489:   } else if (test_mode == 3) {
490:     PetscCall(ex1_3());
491:   } else if (test_mode == 4) {
492:     PetscCall(ex1_4());
493:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown test_mode value, should be 1,2,3,4");
494:   PetscCall(PetscFinalize());
495:   return 0;
496: }

498: /*TEST

500:    build:
501:       requires: !complex double

503:    test:
504:       args: -test_mode 1
505:       filter: grep -v atomic
506:       filter_output: grep -v atomic

508:    test:
509:       suffix: 2
510:       args: -test_mode 2
511:       filter: grep -v atomic
512:       filter_output: grep -v atomic

514:    test:
515:       suffix: 3
516:       args: -test_mode 3
517:       filter: grep -v atomic
518:       filter_output: grep -v atomic
519:       TODO: broken

521:    test:
522:       suffix: 4
523:       args: -test_mode 4
524:       filter: grep -v atomic
525:       filter_output: grep -v atomic

527:    test:
528:       suffix: 5
529:       nsize: 4
530:       args: -test_mode 1
531:       filter: grep -v atomic
532:       filter_output: grep -v atomic

534:    test:
535:       suffix: 6
536:       nsize: 4
537:       args: -test_mode 2
538:       filter: grep -v atomic
539:       filter_output: grep -v atomic

541:    test:
542:       suffix: 7
543:       nsize: 4
544:       args: -test_mode 3
545:       filter: grep -v atomic
546:       filter_output: grep -v atomic
547:       TODO: broken

549:    test:
550:       suffix: 8
551:       nsize: 4
552:       args: -test_mode 4
553:       filter: grep -v atomic
554:       filter_output: grep -v atomic

556: TEST*/