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