Actual source code: swarm_migrate.c
1: #include <petscsf.h>
2: #include <petscdmswarm.h>
3: #include <petscdmda.h>
4: #include <petsc/private/dmswarmimpl.h>
5: #include "../src/dm/impls/swarm/data_bucket.h"
6: #include "../src/dm/impls/swarm/data_ex.h"
8: /*
9: User loads desired location (MPI rank) into field DMSwarm_rank
10: */
11: PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points)
12: {
13: DM_Swarm *swarm = (DM_Swarm*)dm->data;
14: DMSwarmDataEx de;
15: PetscInt p,npoints,*rankval,n_points_recv;
16: PetscMPIInt rank,nrank;
17: void *point_buffer,*recv_points;
18: size_t sizeof_dmswarm_point;
20: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
22: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
23: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
24: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0, &de);
25: DMSwarmDataExTopologyInitialize(de);
26: for (p = 0; p < npoints; ++p) {
27: nrank = rankval[p];
28: if (nrank != rank) {
29: DMSwarmDataExTopologyAddNeighbour(de,nrank);
30: }
31: }
32: DMSwarmDataExTopologyFinalize(de);
33: DMSwarmDataExInitializeSendCount(de);
34: for (p=0; p<npoints; p++) {
35: nrank = rankval[p];
36: if (nrank != rank) {
37: DMSwarmDataExAddToSendCount(de,nrank,1);
38: }
39: }
40: DMSwarmDataExFinalizeSendCount(de);
41: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
42: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
43: for (p=0; p<npoints; p++) {
44: nrank = rankval[p];
45: if (nrank != rank) {
46: /* copy point into buffer */
47: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
48: /* insert point buffer into DMSwarmDataExchanger */
49: DMSwarmDataExPackData(de,nrank,1,point_buffer);
50: }
51: }
52: DMSwarmDataExPackFinalize(de);
53: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
55: if (remove_sent_points) {
56: DMSwarmDataField gfield;
58: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&gfield);
59: DMSwarmDataFieldGetAccess(gfield);
60: DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);
62: /* remove points which left processor */
63: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
64: for (p=0; p<npoints; p++) {
65: nrank = rankval[p];
66: if (nrank != rank) {
67: /* kill point */
68: DMSwarmDataFieldRestoreAccess(gfield);
70: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
72: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
73: DMSwarmDataFieldGetAccess(gfield);
74: DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);
75: p--; /* check replacement point */
76: }
77: }
78: DMSwarmDataFieldRestoreEntries(gfield,(void**)&rankval);
79: DMSwarmDataFieldRestoreAccess(gfield);
80: }
81: DMSwarmDataExBegin(de);
82: DMSwarmDataExEnd(de);
83: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
84: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
85: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
86: for (p=0; p<n_points_recv; p++) {
87: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
89: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
90: }
91: DMSwarmDataExView(de);
92: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
93: DMSwarmDataExDestroy(de);
94: return 0;
95: }
97: PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm,DM dmcell,PetscBool remove_sent_points,PetscInt *npoints_prior_migration)
98: {
99: DM_Swarm *swarm = (DM_Swarm*)dm->data;
100: DMSwarmDataEx de;
101: PetscInt r,p,npoints,*rankval,n_points_recv;
102: PetscMPIInt rank,_rank;
103: const PetscMPIInt *neighbourranks;
104: void *point_buffer,*recv_points;
105: size_t sizeof_dmswarm_point;
106: PetscInt nneighbors;
107: PetscMPIInt mynneigh,*myneigh;
109: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
110: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
111: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
112: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
113: DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);
114: DMSwarmDataExTopologyInitialize(de);
115: for (r=0; r<nneighbors; r++) {
116: _rank = neighbourranks[r];
117: if ((_rank != rank) && (_rank > 0)) {
118: DMSwarmDataExTopologyAddNeighbour(de,_rank);
119: }
120: }
121: DMSwarmDataExTopologyFinalize(de);
122: DMSwarmDataExTopologyGetNeighbours(de,&mynneigh,&myneigh);
123: DMSwarmDataExInitializeSendCount(de);
124: for (p=0; p<npoints; p++) {
125: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
126: for (r=0; r<mynneigh; r++) {
127: _rank = myneigh[r];
128: DMSwarmDataExAddToSendCount(de,_rank,1);
129: }
130: }
131: }
132: DMSwarmDataExFinalizeSendCount(de);
133: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
134: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
135: for (p=0; p<npoints; p++) {
136: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
137: for (r=0; r<mynneigh; r++) {
138: _rank = myneigh[r];
139: /* copy point into buffer */
140: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
141: /* insert point buffer into DMSwarmDataExchanger */
142: DMSwarmDataExPackData(de,_rank,1,point_buffer);
143: }
144: }
145: }
146: DMSwarmDataExPackFinalize(de);
147: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
148: if (remove_sent_points) {
149: DMSwarmDataField PField;
151: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
152: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
153: /* remove points which left processor */
154: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
155: for (p=0; p<npoints; p++) {
156: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
157: /* kill point */
158: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
159: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
160: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
161: p--; /* check replacement point */
162: }
163: }
164: }
165: DMSwarmDataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);
166: DMSwarmDataExBegin(de);
167: DMSwarmDataExEnd(de);
168: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
169: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
170: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
171: for (p=0; p<n_points_recv; p++) {
172: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
174: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
175: }
176: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
177: DMSwarmDataExDestroy(de);
178: return 0;
179: }
181: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm,PetscBool remove_sent_points)
182: {
183: DM_Swarm *swarm = (DM_Swarm*)dm->data;
184: PetscInt p,npoints,npointsg=0,npoints2,npoints2g,*rankval,npoints_prior_migration;
185: PetscSF sfcell = NULL;
186: const PetscSFNode *LA_sfcell;
187: DM dmcell;
188: Vec pos;
189: PetscBool error_check = swarm->migrate_error_on_missing_point;
190: PetscMPIInt size,rank;
192: DMSwarmGetCellDM(dm,&dmcell);
195: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
196: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
198: #if 1
199: {
200: PetscInt *p_cellid;
201: PetscInt npoints_curr,range = 0;
202: PetscSFNode *sf_cells;
204: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
205: PetscMalloc1(npoints_curr, &sf_cells);
207: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
208: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
209: for (p=0; p<npoints_curr; p++) {
211: sf_cells[p].rank = 0;
212: sf_cells[p].index = p_cellid[p];
213: if (p_cellid[p] > range) {
214: range = p_cellid[p];
215: }
216: }
217: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
218: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
220: /* PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell); */
221: PetscSFCreate(PETSC_COMM_SELF,&sfcell);
222: PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
223: }
224: #endif
226: DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
227: DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
228: DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
230: if (error_check) {
231: DMSwarmGetSize(dm,&npointsg);
232: }
233: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
234: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
235: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
236: for (p=0; p<npoints; p++) {
237: rankval[p] = LA_sfcell[p].index;
238: }
239: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
240: PetscSFDestroy(&sfcell);
242: if (size > 1) {
243: DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);
244: } else {
245: DMSwarmDataField PField;
246: PetscInt npoints_curr;
248: /* remove points which the domain */
249: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
250: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
252: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
253: for (p=0; p<npoints_curr; p++) {
254: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
255: /* kill point */
256: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
257: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL); /* you need to update npoints as the list size decreases! */
258: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
259: p--; /* check replacement point */
260: }
261: }
262: DMSwarmGetSize(dm,&npoints_prior_migration);
264: }
266: /* locate points newly received */
267: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
269: #if 0
270: { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
271: PetscScalar *LA_coor;
272: PetscInt bs;
273: DMSwarmDataField PField;
275: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
276: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);
277: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
279: VecDestroy(&pos);
280: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
282: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
283: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
284: for (p=0; p<npoints2; p++) {
285: rankval[p] = LA_sfcell[p].index;
286: }
287: PetscSFDestroy(&sfcell);
288: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
290: /* remove points which left processor */
291: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
292: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
294: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
295: for (p=0; p<npoints2; p++) {
296: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
297: /* kill point */
298: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
299: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
300: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
301: p--; /* check replacement point */
302: }
303: }
304: }
305: #endif
307: { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
308: PetscScalar *LA_coor;
309: PetscInt npoints_from_neighbours,bs;
310: DMSwarmDataField PField;
312: npoints_from_neighbours = npoints2 - npoints_prior_migration;
314: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
315: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints_from_neighbours,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);
317: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
319: VecDestroy(&pos);
320: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
322: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
323: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
324: for (p=0; p<npoints_from_neighbours; p++) {
325: rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
326: }
327: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
328: PetscSFDestroy(&sfcell);
330: /* remove points which left processor */
331: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
332: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
334: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
335: for (p=npoints_prior_migration; p<npoints2; p++) {
336: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
337: /* kill point */
338: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
339: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
340: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
341: p--; /* check replacement point */
342: }
343: }
344: }
346: {
347: PetscInt *p_cellid;
349: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
350: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
351: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
352: for (p=0; p<npoints2; p++) {
353: p_cellid[p] = rankval[p];
354: }
355: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
356: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
357: }
359: /* check for error on removed points */
360: if (error_check) {
361: DMSwarmGetSize(dm,&npoints2g);
363: }
364: return 0;
365: }
367: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points)
368: {
369: return 0;
370: }
372: /*
373: Redundant as this assumes points can only be sent to a single rank
374: */
375: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
376: {
377: DM_Swarm *swarm = (DM_Swarm*)dm->data;
378: DMSwarmDataEx de;
379: PetscInt p,npoints,*rankval,n_points_recv;
380: PetscMPIInt rank,nrank,negrank;
381: void *point_buffer,*recv_points;
382: size_t sizeof_dmswarm_point;
384: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
385: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
386: *globalsize = npoints;
387: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
388: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
389: DMSwarmDataExTopologyInitialize(de);
390: for (p=0; p<npoints; p++) {
391: negrank = rankval[p];
392: if (negrank < 0) {
393: nrank = -negrank - 1;
394: DMSwarmDataExTopologyAddNeighbour(de,nrank);
395: }
396: }
397: DMSwarmDataExTopologyFinalize(de);
398: DMSwarmDataExInitializeSendCount(de);
399: for (p=0; p<npoints; p++) {
400: negrank = rankval[p];
401: if (negrank < 0) {
402: nrank = -negrank - 1;
403: DMSwarmDataExAddToSendCount(de,nrank,1);
404: }
405: }
406: DMSwarmDataExFinalizeSendCount(de);
407: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
408: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
409: for (p=0; p<npoints; p++) {
410: negrank = rankval[p];
411: if (negrank < 0) {
412: nrank = -negrank - 1;
413: rankval[p] = nrank;
414: /* copy point into buffer */
415: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
416: /* insert point buffer into DMSwarmDataExchanger */
417: DMSwarmDataExPackData(de,nrank,1,point_buffer);
418: rankval[p] = negrank;
419: }
420: }
421: DMSwarmDataExPackFinalize(de);
422: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
423: DMSwarmDataExBegin(de);
424: DMSwarmDataExEnd(de);
425: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
426: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
427: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
428: for (p=0; p<n_points_recv; p++) {
429: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
431: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
432: }
433: DMSwarmDataExView(de);
434: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
435: DMSwarmDataExDestroy(de);
436: return 0;
437: }
439: typedef struct {
440: PetscMPIInt owner_rank;
441: PetscReal min[3],max[3];
442: } CollectBBox;
444: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
445: {
446: DM_Swarm * swarm = (DM_Swarm*)dm->data;
447: DMSwarmDataEx de;
448: PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
449: PetscMPIInt rank,nrank;
450: void *point_buffer,*recv_points;
451: size_t sizeof_dmswarm_point,sizeof_bbox_ctx;
452: PetscBool isdmda;
453: CollectBBox *bbox,*recv_bbox;
454: const PetscMPIInt *dmneighborranks;
455: DM dmcell;
457: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
459: DMSwarmGetCellDM(dm,&dmcell);
461: isdmda = PETSC_FALSE;
462: PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);
465: DMGetDimension(dm,&dim);
466: sizeof_bbox_ctx = sizeof(CollectBBox);
467: PetscMalloc1(1,&bbox);
468: bbox->owner_rank = rank;
470: /* compute the bounding box based on the overlapping / stenctil size */
471: {
472: Vec lcoor;
474: DMGetCoordinatesLocal(dmcell,&lcoor);
475: if (dim >= 1) {
476: VecStrideMin(lcoor,0,NULL,&bbox->min[0]);
477: VecStrideMax(lcoor,0,NULL,&bbox->max[0]);
478: }
479: if (dim >= 2) {
480: VecStrideMin(lcoor,1,NULL,&bbox->min[1]);
481: VecStrideMax(lcoor,1,NULL,&bbox->max[1]);
482: }
483: if (dim == 3) {
484: VecStrideMin(lcoor,2,NULL,&bbox->min[2]);
485: VecStrideMax(lcoor,2,NULL,&bbox->max[2]);
486: }
487: }
488: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
489: *globalsize = npoints;
490: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
491: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
492: /* use DMDA neighbours */
493: DMDAGetNeighbors(dmcell,&dmneighborranks);
494: if (dim == 1) {
495: neighbour_cells = 3;
496: } else if (dim == 2) {
497: neighbour_cells = 9;
498: } else {
499: neighbour_cells = 27;
500: }
501: DMSwarmDataExTopologyInitialize(de);
502: for (p=0; p<neighbour_cells; p++) {
503: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) {
504: DMSwarmDataExTopologyAddNeighbour(de,dmneighborranks[p]);
505: }
506: }
507: DMSwarmDataExTopologyFinalize(de);
508: DMSwarmDataExInitializeSendCount(de);
509: for (p=0; p<neighbour_cells; p++) {
510: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) {
511: DMSwarmDataExAddToSendCount(de,dmneighborranks[p],1);
512: }
513: }
514: DMSwarmDataExFinalizeSendCount(de);
515: /* send bounding boxes */
516: DMSwarmDataExPackInitialize(de,sizeof_bbox_ctx);
517: for (p=0; p<neighbour_cells; p++) {
518: nrank = dmneighborranks[p];
519: if ((nrank >= 0) && (nrank != rank)) {
520: /* insert bbox buffer into DMSwarmDataExchanger */
521: DMSwarmDataExPackData(de,nrank,1,bbox);
522: }
523: }
524: DMSwarmDataExPackFinalize(de);
525: /* recv bounding boxes */
526: DMSwarmDataExBegin(de);
527: DMSwarmDataExEnd(de);
528: DMSwarmDataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);
529: /* Wrong, should not be using PETSC_COMM_WORLD */
530: for (p=0; p<n_bbox_recv; p++) {
531: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank,
532: (double)recv_bbox[p].min[0],(double)recv_bbox[p].max[0],(double)recv_bbox[p].min[1],(double)recv_bbox[p].max[1]));
533: }
534: PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
535: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
536: DMSwarmDataExInitializeSendCount(de);
537: for (pk=0; pk<n_bbox_recv; pk++) {
538: PetscReal *array_x,*array_y;
540: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
541: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
542: for (p=0; p<npoints; p++) {
543: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
544: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
545: DMSwarmDataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);
546: }
547: }
548: }
549: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
550: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
551: }
552: DMSwarmDataExFinalizeSendCount(de);
553: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
554: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
555: for (pk=0; pk<n_bbox_recv; pk++) {
556: PetscReal *array_x,*array_y;
558: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
559: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
560: for (p=0; p<npoints; p++) {
561: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
562: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
563: /* copy point into buffer */
564: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
565: /* insert point buffer into DMSwarmDataExchanger */
566: DMSwarmDataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);
567: }
568: }
569: }
570: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
571: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
572: }
573: DMSwarmDataExPackFinalize(de);
574: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
575: DMSwarmDataExBegin(de);
576: DMSwarmDataExEnd(de);
577: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
578: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
579: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
580: for (p=0; p<n_points_recv; p++) {
581: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
583: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
584: }
585: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
586: PetscFree(bbox);
587: DMSwarmDataExView(de);
588: DMSwarmDataExDestroy(de);
589: return 0;
590: }
592: /* General collection when no order, or neighbour information is provided */
593: /*
594: User provides context and collect() method
595: Broadcast user context
597: for each context / rank {
598: collect(swarm,context,n,list)
599: }
600: */
601: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
602: {
603: DM_Swarm *swarm = (DM_Swarm*)dm->data;
604: DMSwarmDataEx de;
605: PetscInt p,r,npoints,n_points_recv;
606: PetscMPIInt size,rank;
607: void *point_buffer,*recv_points;
608: void *ctxlist;
609: PetscInt *n2collect,**collectlist;
610: size_t sizeof_dmswarm_point;
612: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
613: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
614: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
615: *globalsize = npoints;
616: /* Broadcast user context */
617: PetscMalloc(ctx_size*size,&ctxlist);
618: MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));
619: PetscMalloc1(size,&n2collect);
620: PetscMalloc1(size,&collectlist);
621: for (r=0; r<size; r++) {
622: PetscInt _n2collect;
623: PetscInt *_collectlist;
624: void *_ctx_r;
626: _n2collect = 0;
627: _collectlist = NULL;
628: if (r != rank) { /* don't collect data from yourself */
629: _ctx_r = (void*)( (char*)ctxlist + r * ctx_size);
630: collect(dm,_ctx_r,&_n2collect,&_collectlist);
631: }
632: n2collect[r] = _n2collect;
633: collectlist[r] = _collectlist;
634: }
635: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
636: /* Define topology */
637: DMSwarmDataExTopologyInitialize(de);
638: for (r=0; r<size; r++) {
639: if (n2collect[r] > 0) {
640: DMSwarmDataExTopologyAddNeighbour(de,(PetscMPIInt)r);
641: }
642: }
643: DMSwarmDataExTopologyFinalize(de);
644: /* Define send counts */
645: DMSwarmDataExInitializeSendCount(de);
646: for (r=0; r<size; r++) {
647: if (n2collect[r] > 0) {
648: DMSwarmDataExAddToSendCount(de,r,n2collect[r]);
649: }
650: }
651: DMSwarmDataExFinalizeSendCount(de);
652: /* Pack data */
653: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
654: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
655: for (r=0; r<size; r++) {
656: for (p=0; p<n2collect[r]; p++) {
657: DMSwarmDataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);
658: /* insert point buffer into the data exchanger */
659: DMSwarmDataExPackData(de,r,1,point_buffer);
660: }
661: }
662: DMSwarmDataExPackFinalize(de);
663: /* Scatter */
664: DMSwarmDataExBegin(de);
665: DMSwarmDataExEnd(de);
666: /* Collect data in DMSwarm container */
667: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
668: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
669: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
670: for (p=0; p<n_points_recv; p++) {
671: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
673: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
674: }
675: /* Release memory */
676: for (r=0; r<size; r++) {
677: if (collectlist[r]) PetscFree(collectlist[r]);
678: }
679: PetscFree(collectlist);
680: PetscFree(n2collect);
681: PetscFree(ctxlist);
682: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
683: DMSwarmDataExView(de);
684: DMSwarmDataExDestroy(de);
685: return 0;
686: }