Actual source code: swarm_ex1.c


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

181:   PetscFunctionBegin;
182:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
183:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
184:   overlap = 2;
185:   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));
186:   PetscCall(DMSetFromOptions(dmcell));
187:   PetscCall(DMSetUp(dmcell));
188:   PetscCall(DMDASetUniformCoordinates(dmcell, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
189:   PetscCall(DMDAGetCorners(dmcell, &is, &js, NULL, &ni, &nj, NULL));
190:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
191:   PetscCall(DMSetType(dms, DMSWARM));
192:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));
193:   PetscCall(DMSwarmSetCellDM(dms, dmcell));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

499: /*TEST

501:    build:
502:       requires: !complex double

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

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

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

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

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

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

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

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

557: TEST*/