xref: /petsc/src/dm/tutorials/swarm_ex1.c (revision 3307d110e72ee4e6d2468971620073eb5ff93529)
1 
2 static char help[] = "Tests DMSwarm\n\n";
3 
4 #include <petscdm.h>
5 #include <petscdmda.h>
6 #include <petscdmswarm.h>
7 
8 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM,PetscErrorCode (*)(DM,void*,PetscInt*,PetscInt**),size_t,void*,PetscInt*);
9 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM,PetscInt*);
10 
11 PetscErrorCode ex1_1(void)
12 {
13   DM             dms;
14   Vec            x;
15   PetscMPIInt    rank,size;
16   PetscInt       p;
17 
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 wuth 4 MPI ranks");
22 
23   PetscCall(DMCreate(PETSC_COMM_WORLD,&dms));
24   PetscCall(DMSetType(dms,DMSWARM));
25   PetscCall(PetscObjectSetName((PetscObject) dms, "Particles"));
26 
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));
33 
34   {
35     PetscReal *array;
36     PetscCall(DMSwarmGetField(dms,"viscosity",NULL,NULL,(void**)&array));
37     for (p=0; p<5+rank; p++) {
38       array[p] = 11.1 + p*0.1 + rank*100.0;
39     }
40     PetscCall(DMSwarmRestoreField(dms,"viscosity",NULL,NULL,(void**)&array));
41   }
42 
43   {
44     PetscReal *array;
45     PetscCall(DMSwarmGetField(dms,"strain",NULL,NULL,(void**)&array));
46     for (p=0; p<5+rank; p++) {
47       array[p] = 2.0e-2 + p*0.001 + rank*1.0;
48     }
49     PetscCall(DMSwarmRestoreField(dms,"strain",NULL,NULL,(void**)&array));
50   }
51 
52   PetscCall(DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x));
53   PetscCall(DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x));
54 
55   PetscCall(DMSwarmVectorDefineField(dms,"strain"));
56   PetscCall(DMCreateGlobalVector(dms,&x));
57   PetscCall(VecDestroy(&x));
58 
59   {
60     PetscInt    *rankval;
61     PetscInt    npoints[2],npoints_orig[2];
62 
63     PetscCall(DMSwarmGetLocalSize(dms,&npoints_orig[0]));
64     PetscCall(DMSwarmGetSize(dms,&npoints_orig[1]));
65     PetscCall(DMSwarmGetField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval));
66     if ((rank == 0) && (size > 1)) {
67       rankval[0] = 1;
68       rankval[3] = 1;
69     }
70     if (rank == 3) {
71       rankval[2] = 1;
72     }
73     PetscCall(DMSwarmRestoreField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval));
74     PetscCall(DMSwarmMigrate(dms,PETSC_TRUE));
75     PetscCall(DMSwarmGetLocalSize(dms,&npoints[0]));
76     PetscCall(DMSwarmGetSize(dms,&npoints[1]));
77     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
78     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]));
79     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
80     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
81   }
82   {
83     PetscCall(DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x));
84     PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
85     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x));
86   }
87   {
88     PetscCall(DMSwarmCreateGlobalVectorFromField(dms,"strain",&x));
89     PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
90     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms,"strain",&x));
91   }
92 
93   PetscCall(DMDestroy(&dms));
94   PetscFunctionReturn(0);
95 }
96 
97 PetscErrorCode ex1_2(void)
98 {
99   DM             dms;
100   Vec            x;
101   PetscMPIInt    rank,size;
102   PetscInt       p;
103 
104   PetscFunctionBegin;
105   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
106   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
107   PetscCall(DMCreate(PETSC_COMM_WORLD,&dms));
108   PetscCall(DMSetType(dms,DMSWARM));
109   PetscCall(PetscObjectSetName((PetscObject) dms, "Particles"));
110   PetscCall(DMSwarmInitializeFieldRegister(dms));
111 
112   PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL));
113   PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"strain",1,PETSC_REAL));
114   PetscCall(DMSwarmFinalizeFieldRegister(dms));
115   PetscCall(DMSwarmSetLocalSizes(dms,5+rank,4));
116   PetscCall(DMView(dms,PETSC_VIEWER_STDOUT_WORLD));
117   {
118     PetscReal *array;
119     PetscCall(DMSwarmGetField(dms,"viscosity",NULL,NULL,(void**)&array));
120     for (p=0; p<5+rank; p++) {
121       array[p] = 11.1 + p*0.1 + rank*100.0;
122     }
123     PetscCall(DMSwarmRestoreField(dms,"viscosity",NULL,NULL,(void**)&array));
124   }
125   {
126     PetscReal *array;
127     PetscCall(DMSwarmGetField(dms,"strain",NULL,NULL,(void**)&array));
128     for (p=0; p<5+rank; p++) {
129       array[p] = 2.0e-2 + p*0.001 + rank*1.0;
130     }
131     PetscCall(DMSwarmRestoreField(dms,"strain",NULL,NULL,(void**)&array));
132   }
133   {
134     PetscInt    *rankval;
135     PetscInt    npoints[2],npoints_orig[2];
136 
137     PetscCall(DMSwarmGetLocalSize(dms,&npoints_orig[0]));
138     PetscCall(DMSwarmGetSize(dms,&npoints_orig[1]));
139     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
140     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n",rank,npoints_orig[0],npoints_orig[1]));
141     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
142     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
143 
144     PetscCall(DMSwarmGetField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval));
145 
146     if (rank == 1) {
147       rankval[0] = -1;
148     }
149     if (rank == 2) {
150       rankval[1] = -1;
151     }
152     if (rank == 3) {
153       rankval[3] = -1;
154       rankval[4] = -1;
155     }
156     PetscCall(DMSwarmRestoreField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval));
157     PetscCall(DMSwarmCollectViewCreate(dms));
158     PetscCall(DMSwarmGetLocalSize(dms,&npoints[0]));
159     PetscCall(DMSwarmGetSize(dms,&npoints[1]));
160     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
161     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n",rank,npoints[0],npoints[1]));
162     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
163     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
164 
165     PetscCall(DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x));
166     PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
167     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x));
168 
169     PetscCall(DMSwarmCollectViewDestroy(dms));
170     PetscCall(DMSwarmGetLocalSize(dms,&npoints[0]));
171     PetscCall(DMSwarmGetSize(dms,&npoints[1]));
172     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
173     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after_v(%" PetscInt_FMT ",%" PetscInt_FMT ")\n",rank,npoints[0],npoints[1]));
174     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
175     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
176 
177     PetscCall(DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x));
178     PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
179     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x));
180   }
181   PetscCall(DMDestroy(&dms));
182   PetscFunctionReturn(0);
183 }
184 
185 /*
186  splot "c-rank0.gp","c-rank1.gp","c-rank2.gp","c-rank3.gp"
187 */
188 PetscErrorCode ex1_3(void)
189 {
190   DM             dms;
191   PetscMPIInt    rank,size;
192   PetscInt       is,js,ni,nj,overlap;
193   DM             dmcell;
194 
195   PetscFunctionBegin;
196   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
197   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
198   overlap = 2;
199   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));
200   PetscCall(DMSetFromOptions(dmcell));
201   PetscCall(DMSetUp(dmcell));
202   PetscCall(DMDASetUniformCoordinates(dmcell,-1.0,1.0,-1.0,1.0,-1.0,1.0));
203   PetscCall(DMDAGetCorners(dmcell,&is,&js,NULL,&ni,&nj,NULL));
204   PetscCall(DMCreate(PETSC_COMM_WORLD,&dms));
205   PetscCall(DMSetType(dms,DMSWARM));
206   PetscCall(PetscObjectSetName((PetscObject) dms, "Particles"));
207   PetscCall(DMSwarmSetCellDM(dms,dmcell));
208 
209   /* load in data types */
210   PetscCall(DMSwarmInitializeFieldRegister(dms));
211   PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL));
212   PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"coorx",1,PETSC_REAL));
213   PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"coory",1,PETSC_REAL));
214   PetscCall(DMSwarmFinalizeFieldRegister(dms));
215   PetscCall(DMSwarmSetLocalSizes(dms,ni*nj*4,4));
216   PetscCall(DMView(dms,PETSC_VIEWER_STDOUT_WORLD));
217 
218   /* set values within the swarm */
219   {
220     PetscReal  *array_x,*array_y;
221     PetscInt   npoints,i,j,cnt;
222     DMDACoor2d **LA_coor;
223     Vec        coor;
224     DM         dmcellcdm;
225 
226     PetscCall(DMGetCoordinateDM(dmcell,&dmcellcdm));
227     PetscCall(DMGetCoordinates(dmcell,&coor));
228     PetscCall(DMDAVecGetArray(dmcellcdm,coor,&LA_coor));
229     PetscCall(DMSwarmGetLocalSize(dms,&npoints));
230     PetscCall(DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x));
231     PetscCall(DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y));
232     cnt = 0;
233     for (j=js; j<js+nj; j++) {
234       for (i=is; i<is+ni; i++) {
235         PetscReal xp,yp;
236         xp = PetscRealPart(LA_coor[j][i].x);
237         yp = PetscRealPart(LA_coor[j][i].y);
238         array_x[4*cnt+0] = xp - 0.05; if (array_x[4*cnt+0] < -1.0) { array_x[4*cnt+0] = -1.0+1.0e-12; }
239         array_x[4*cnt+1] = xp + 0.05; if (array_x[4*cnt+1] > 1.0)  { array_x[4*cnt+1] =  1.0-1.0e-12; }
240         array_x[4*cnt+2] = xp - 0.05; if (array_x[4*cnt+2] < -1.0) { array_x[4*cnt+2] = -1.0+1.0e-12; }
241         array_x[4*cnt+3] = xp + 0.05; if (array_x[4*cnt+3] > 1.0)  { array_x[4*cnt+3] =  1.0-1.0e-12; }
242 
243         array_y[4*cnt+0] = yp - 0.05; if (array_y[4*cnt+0] < -1.0) { array_y[4*cnt+0] = -1.0+1.0e-12; }
244         array_y[4*cnt+1] = yp - 0.05; if (array_y[4*cnt+1] < -1.0) { array_y[4*cnt+1] = -1.0+1.0e-12; }
245         array_y[4*cnt+2] = yp + 0.05; if (array_y[4*cnt+2] > 1.0)  { array_y[4*cnt+2] =  1.0-1.0e-12; }
246         array_y[4*cnt+3] = yp + 0.05; if (array_y[4*cnt+3] > 1.0)  { array_y[4*cnt+3] =  1.0-1.0e-12; }
247         cnt++;
248       }
249     }
250     PetscCall(DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y));
251     PetscCall(DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x));
252     PetscCall(DMDAVecRestoreArray(dmcellcdm,coor,&LA_coor));
253   }
254   {
255     PetscInt    npoints[2],npoints_orig[2],ng;
256 
257     PetscCall(DMSwarmGetLocalSize(dms,&npoints_orig[0]));
258     PetscCall(DMSwarmGetSize(dms,&npoints_orig[1]));
259     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
260     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n",rank,npoints_orig[0],npoints_orig[1]));
261     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
262     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
263     PetscCall(DMSwarmCollect_DMDABoundingBox(dms,&ng));
264 
265     PetscCall(DMSwarmGetLocalSize(dms,&npoints[0]));
266     PetscCall(DMSwarmGetSize(dms,&npoints[1]));
267     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
268     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n",rank,npoints[0],npoints[1]));
269     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
270     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
271   }
272   {
273     PetscReal *array_x,*array_y;
274     PetscInt  npoints,p;
275     FILE      *fp = NULL;
276     char      name[PETSC_MAX_PATH_LEN];
277 
278     PetscCall(PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"c-rank%d.gp",rank));
279     fp = fopen(name,"w");
280     PetscCheck(fp,PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file %s",name);
281     PetscCall(DMSwarmGetLocalSize(dms,&npoints));
282     PetscCall(DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x));
283     PetscCall(DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y));
284     for (p=0; p<npoints; p++) {
285       fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array_x[p],array_y[p],(double)rank);
286     }
287     PetscCall(DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y));
288     PetscCall(DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x));
289     fclose(fp);
290   }
291   PetscCall(DMDestroy(&dmcell));
292   PetscCall(DMDestroy(&dms));
293   PetscFunctionReturn(0);
294 }
295 
296 typedef struct {
297   PetscReal cx[2];
298   PetscReal radius;
299 } CollectZoneCtx;
300 
301 PetscErrorCode collect_zone(DM dm,void *ctx,PetscInt *nfound,PetscInt **foundlist)
302 {
303   CollectZoneCtx *zone = (CollectZoneCtx*)ctx;
304   PetscInt       p,npoints;
305   PetscReal      *array_x,*array_y,r2;
306   PetscInt       p2collect,*plist;
307   PetscMPIInt    rank;
308 
309   PetscFunctionBegin;
310   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
311   PetscCall(DMSwarmGetLocalSize(dm,&npoints));
312   PetscCall(DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x));
313   PetscCall(DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y));
314 
315   r2 = zone->radius * zone->radius;
316   p2collect = 0;
317   for (p=0; p<npoints; p++) {
318     PetscReal sep2;
319 
320     sep2  = (array_x[p] - zone->cx[0])*(array_x[p] - zone->cx[0]);
321     sep2 += (array_y[p] - zone->cx[1])*(array_y[p] - zone->cx[1]);
322     if (sep2 < r2) {
323       p2collect++;
324     }
325   }
326 
327   PetscCall(PetscMalloc1(p2collect+1,&plist));
328   p2collect = 0;
329   for (p=0; p<npoints; p++) {
330     PetscReal sep2;
331 
332     sep2  = (array_x[p] - zone->cx[0])*(array_x[p] - zone->cx[0]);
333     sep2 += (array_y[p] - zone->cx[1])*(array_y[p] - zone->cx[1]);
334     if (sep2 < r2) {
335       plist[p2collect] = p;
336       p2collect++;
337     }
338   }
339   PetscCall(DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y));
340   PetscCall(DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x));
341 
342   *nfound = p2collect;
343   *foundlist = plist;
344   PetscFunctionReturn(0);
345 }
346 
347 PetscErrorCode ex1_4(void)
348 {
349   DM             dms;
350   PetscMPIInt    rank,size;
351   PetscInt       is,js,ni,nj,overlap,nn;
352   DM             dmcell;
353   CollectZoneCtx *zone;
354   PetscReal      dx;
355 
356   PetscFunctionBegin;
357   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
358   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
359   nn = 101;
360   dx = 2.0/ (PetscReal)(nn-1);
361   overlap = 0;
362   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));
363   PetscCall(DMSetFromOptions(dmcell));
364   PetscCall(DMSetUp(dmcell));
365   PetscCall(DMDASetUniformCoordinates(dmcell,-1.0,1.0,-1.0,1.0,-1.0,1.0));
366   PetscCall(DMDAGetCorners(dmcell,&is,&js,NULL,&ni,&nj,NULL));
367   PetscCall(DMCreate(PETSC_COMM_WORLD,&dms));
368   PetscCall(DMSetType(dms,DMSWARM));
369   PetscCall(PetscObjectSetName((PetscObject) dms, "Particles"));
370 
371   /* load in data types */
372   PetscCall(DMSwarmInitializeFieldRegister(dms));
373   PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL));
374   PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"coorx",1,PETSC_REAL));
375   PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"coory",1,PETSC_REAL));
376   PetscCall(DMSwarmFinalizeFieldRegister(dms));
377   PetscCall(DMSwarmSetLocalSizes(dms,ni*nj*4,4));
378   PetscCall(DMView(dms,PETSC_VIEWER_STDOUT_WORLD));
379 
380   /* set values within the swarm */
381   {
382     PetscReal  *array_x,*array_y;
383     PetscInt   npoints,i,j,cnt;
384     DMDACoor2d **LA_coor;
385     Vec        coor;
386     DM         dmcellcdm;
387 
388     PetscCall(DMGetCoordinateDM(dmcell,&dmcellcdm));
389     PetscCall(DMGetCoordinates(dmcell,&coor));
390     PetscCall(DMDAVecGetArray(dmcellcdm,coor,&LA_coor));
391     PetscCall(DMSwarmGetLocalSize(dms,&npoints));
392     PetscCall(DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x));
393     PetscCall(DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y));
394     cnt = 0;
395     for (j=js; j<js+nj; j++) {
396       for (i=is; i<is+ni; i++) {
397         PetscReal xp,yp;
398 
399         xp = PetscRealPart(LA_coor[j][i].x);
400         yp = PetscRealPart(LA_coor[j][i].y);
401         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; }*/
402         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; }*/
403         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; }*/
404         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; }*/
405         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; }*/
406         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; }*/
407         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; }*/
408         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; }*/
409         cnt++;
410       }
411     }
412     PetscCall(DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y));
413     PetscCall(DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x));
414     PetscCall(DMDAVecRestoreArray(dmcellcdm,coor,&LA_coor));
415   }
416   PetscCall(PetscMalloc1(1,&zone));
417   if (size == 4) {
418     if (rank == 0) {
419       zone->cx[0] = 0.5;
420       zone->cx[1] = 0.5;
421       zone->radius = 0.3;
422     }
423     if (rank == 1) {
424       zone->cx[0] = -0.5;
425       zone->cx[1] = 0.5;
426       zone->radius = 0.25;
427     }
428     if (rank == 2) {
429       zone->cx[0] = 0.5;
430       zone->cx[1] = -0.5;
431       zone->radius = 0.2;
432     }
433     if (rank == 3) {
434       zone->cx[0] = -0.5;
435       zone->cx[1] = -0.5;
436       zone->radius = 0.1;
437     }
438   } else {
439     if (rank == 0) {
440       zone->cx[0] = 0.5;
441       zone->cx[1] = 0.5;
442       zone->radius = 0.8;
443     } else {
444       zone->cx[0] = 10.0;
445       zone->cx[1] = 10.0;
446       zone->radius = 0.0;
447     }
448   }
449   {
450     PetscInt    npoints[2],npoints_orig[2],ng;
451 
452     PetscCall(DMSwarmGetLocalSize(dms,&npoints_orig[0]));
453     PetscCall(DMSwarmGetSize(dms,&npoints_orig[1]));
454     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
455     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n",rank,npoints_orig[0],npoints_orig[1]));
456     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
457     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
458     PetscCall(DMSwarmCollect_General(dms,collect_zone,sizeof(CollectZoneCtx),zone,&ng));
459     PetscCall(DMSwarmGetLocalSize(dms,&npoints[0]));
460     PetscCall(DMSwarmGetSize(dms,&npoints[1]));
461     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
462     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n",rank,npoints[0],npoints[1]));
463     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
464     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
465   }
466   {
467     PetscReal *array_x,*array_y;
468     PetscInt  npoints,p;
469     FILE      *fp = NULL;
470     char      name[PETSC_MAX_PATH_LEN];
471 
472     PetscCall(PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"c-rank%d.gp",rank));
473     fp = fopen(name,"w");
474     PetscCheck(fp,PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file %s",name);
475     PetscCall(DMSwarmGetLocalSize(dms,&npoints));
476     PetscCall(DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x));
477     PetscCall(DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y));
478     for (p=0; p<npoints; p++) {
479       fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array_x[p],array_y[p],(double)rank);
480     }
481     PetscCall(DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y));
482     PetscCall(DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x));
483     fclose(fp);
484   }
485   PetscCall(DMDestroy(&dmcell));
486   PetscCall(DMDestroy(&dms));
487   PetscCall(PetscFree(zone));
488   PetscFunctionReturn(0);
489 }
490 
491 int main(int argc,char **argv)
492 {
493   PetscInt       test_mode = 4;
494 
495   PetscCall(PetscInitialize(&argc,&argv,(char*)0,help));
496   PetscCall(PetscOptionsGetInt(NULL,NULL,"-test_mode",&test_mode,NULL));
497   if (test_mode == 1) {
498     PetscCall(ex1_1());
499   } else if (test_mode == 2) {
500     PetscCall(ex1_2());
501   } else if (test_mode == 3) {
502     PetscCall(ex1_3());
503   } else if (test_mode == 4) {
504     PetscCall(ex1_4());
505   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown test_mode value, should be 1,2,3,4");
506   PetscCall(PetscFinalize());
507   return 0;
508 }
509 
510 /*TEST
511 
512    build:
513       requires: !complex double
514 
515    test:
516       args: -test_mode 1
517       filter: grep -v atomic
518       filter_output: grep -v atomic
519 
520    test:
521       suffix: 2
522       args: -test_mode 2
523       filter: grep -v atomic
524       filter_output: grep -v atomic
525 
526    test:
527       suffix: 3
528       args: -test_mode 3
529       filter: grep -v atomic
530       filter_output: grep -v atomic
531       TODO: broken
532 
533    test:
534       suffix: 4
535       args: -test_mode 4
536       filter: grep -v atomic
537       filter_output: grep -v atomic
538 
539    test:
540       suffix: 5
541       nsize: 4
542       args: -test_mode 1
543       filter: grep -v atomic
544       filter_output: grep -v atomic
545 
546    test:
547       suffix: 6
548       nsize: 4
549       args: -test_mode 2
550       filter: grep -v atomic
551       filter_output: grep -v atomic
552 
553    test:
554       suffix: 7
555       nsize: 4
556       args: -test_mode 3
557       filter: grep -v atomic
558       filter_output: grep -v atomic
559       TODO: broken
560 
561    test:
562       suffix: 8
563       nsize: 4
564       args: -test_mode 4
565       filter: grep -v atomic
566       filter_output: grep -v atomic
567 
568 TEST*/
569