Actual source code: swarm_migrate.c
petsc-3.8.4 2018-03-24
1: #include <petscsf.h>
2: #include <petscdmswarm.h>
3: #include <petscdmda.h>
4: #include <petsc/private/dmswarmimpl.h>
5: #include "data_bucket.h"
6: #include "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;
15: DataEx de;
16: PetscInt p,npoints,*rankval,n_points_recv;
17: PetscMPIInt rank,nrank;
18: void *point_buffer,*recv_points;
19: size_t sizeof_dmswarm_point;
22: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
24: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
25: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
26: DataExCreate(PetscObjectComm((PetscObject)dm),0, &de);
27: DataExTopologyInitialize(de);
28: for (p = 0; p < npoints; ++p) {
29: nrank = rankval[p];
30: if (nrank != rank) {
31: DataExTopologyAddNeighbour(de,nrank);
32: }
33: }
34: DataExTopologyFinalize(de);
35: DataExInitializeSendCount(de);
36: for (p=0; p<npoints; p++) {
37: nrank = rankval[p];
38: if (nrank != rank) {
39: DataExAddToSendCount(de,nrank,1);
40: }
41: }
42: DataExFinalizeSendCount(de);
43: DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
44: DataExPackInitialize(de,sizeof_dmswarm_point);
45: for (p=0; p<npoints; p++) {
46: nrank = rankval[p];
47: if (nrank != rank) {
48: /* copy point into buffer */
49: DataBucketFillPackedArray(swarm->db,p,point_buffer);
50: /* insert point buffer into DataExchanger */
51: DataExPackData(de,nrank,1,point_buffer);
52: }
53: }
54: DataExPackFinalize(de);
55: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
57: if (remove_sent_points) {
58: DataField gfield;
60: DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&gfield);
61: DataFieldGetAccess(gfield);
62: DataFieldGetEntries(gfield,(void**)&rankval);
64: /* remove points which left processor */
65: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
66: for (p=0; p<npoints; p++) {
67: nrank = rankval[p];
68: if (nrank != rank) {
69: /* kill point */
70: DataFieldRestoreAccess(gfield);
71:
72: DataBucketRemovePointAtIndex(swarm->db,p);
74: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
75: DataFieldGetAccess(gfield);
76: DataFieldGetEntries(gfield,(void**)&rankval);
77: p--; /* check replacement point */
78: }
79: }
80: DataFieldRestoreEntries(gfield,(void**)&rankval);
81: DataFieldRestoreAccess(gfield);
82: }
83: DataExBegin(de);
84: DataExEnd(de);
85: DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
86: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
87: DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);
88: for (p=0; p<n_points_recv; p++) {
89: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
91: DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
92: }
93: DataExView(de);
94: DataBucketDestroyPackedArray(swarm->db,&point_buffer);
95: DataExDestroy(de);
96: return(0);
97: }
99: PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm,DM dmcell,PetscBool remove_sent_points,PetscInt *npoints_prior_migration)
100: {
101: DM_Swarm *swarm = (DM_Swarm*)dm->data;
103: DataEx de;
104: PetscInt r,p,npoints,*rankval,n_points_recv;
105: PetscMPIInt rank,_rank;
106: const PetscMPIInt *neighbourranks;
107: void *point_buffer,*recv_points;
108: size_t sizeof_dmswarm_point;
109: PetscInt nneighbors;
110: PetscMPIInt mynneigh,*myneigh;
113: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
114: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
115: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
116: DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
117: DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);
118: DataExTopologyInitialize(de);
119: for (r=0; r<nneighbors; r++) {
120: _rank = neighbourranks[r];
121: if ((_rank != rank) && (_rank > 0)) {
122: DataExTopologyAddNeighbour(de,_rank);
123: }
124: }
125: DataExTopologyFinalize(de);
126: DataExTopologyGetNeighbours(de,&mynneigh,&myneigh);
127: DataExInitializeSendCount(de);
128: for (p=0; p<npoints; p++) {
129: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
130: for (r=0; r<mynneigh; r++) {
131: _rank = myneigh[r];
132: DataExAddToSendCount(de,_rank,1);
133: }
134: }
135: }
136: DataExFinalizeSendCount(de);
137: DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
138: DataExPackInitialize(de,sizeof_dmswarm_point);
139: for (p=0; p<npoints; p++) {
140: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
141: for (r=0; r<mynneigh; r++) {
142: _rank = myneigh[r];
143: /* copy point into buffer */
144: DataBucketFillPackedArray(swarm->db,p,point_buffer);
145: /* insert point buffer into DataExchanger */
146: DataExPackData(de,_rank,1,point_buffer);
147: }
148: }
149: }
150: DataExPackFinalize(de);
151: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
152: if (remove_sent_points) {
153: DataField PField;
155: DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
156: DataFieldGetEntries(PField,(void**)&rankval);
157: /* remove points which left processor */
158: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
159: for (p=0; p<npoints; p++) {
160: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
161: /* kill point */
162: DataBucketRemovePointAtIndex(swarm->db,p);
163: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
164: DataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
165: p--; /* check replacement point */
166: }
167: }
168: }
169: DataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);
170: DataExBegin(de);
171: DataExEnd(de);
172: DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
173: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
174: DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);
175: for (p=0; p<n_points_recv; p++) {
176: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
178: DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
179: }
180: DataBucketDestroyPackedArray(swarm->db,&point_buffer);
181: DataExDestroy(de);
182: return(0);
183: }
185: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm,PetscBool remove_sent_points)
186: {
187: DM_Swarm *swarm = (DM_Swarm*)dm->data;
189: PetscInt p,npoints,npointsg=0,npoints2,npoints2g,*rankval,npoints_prior_migration;
190: PetscSF sfcell = NULL;
191: const PetscSFNode *LA_sfcell;
192: DM dmcell;
193: Vec pos;
194: PetscBool error_check = swarm->migrate_error_on_missing_point;
195: PetscMPIInt commsize,rank;
198: DMSwarmGetCellDM(dm,&dmcell);
199: if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
201: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);
202: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
204: #if 1
205: {
206: PetscInt *p_cellid;
207: PetscInt npoints_curr,range = 0;
208: PetscSFNode *sf_cells;
210:
211: DataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
212: PetscMalloc1(npoints_curr, &sf_cells);
214: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
215: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
216: for (p=0; p<npoints_curr; p++) {
218: sf_cells[p].rank = 0;
219: sf_cells[p].index = p_cellid[p];
220: if (p_cellid[p] > range) {
221: range = p_cellid[p];
222: }
223: }
224: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
225: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
227: /*PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell);*/
228: PetscSFCreate(PETSC_COMM_SELF,&sfcell);
229: PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
230: }
231: #endif
232:
233: DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
234: DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
235: DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
237: if (error_check) {
238: DMSwarmGetSize(dm,&npointsg);
239: }
240: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
241: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
242: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
243: for (p=0; p<npoints; p++) {
244: rankval[p] = LA_sfcell[p].index;
245: }
246: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
247: PetscSFDestroy(&sfcell);
249: if (commsize > 1) {
250: DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);
251: } else {
252: DataField PField;
253: PetscInt npoints_curr;
254:
255: /* remove points which the domain */
256: DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
257: DataFieldGetEntries(PField,(void**)&rankval);
258:
259: DataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
260: for (p=0; p<npoints_curr; p++) {
261: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
262: /* kill point */
263: DataBucketRemovePointAtIndex(swarm->db,p);
264: DataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL); /* you need to update npoints as the list size decreases! */
265: DataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
266: p--; /* check replacement point */
267: }
268: }
269: DMSwarmGetSize(dm,&npoints_prior_migration);
270:
271: }
273: /* locate points newly recevied */
274: DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
275:
276: #if 0
277: { /* safe alternative - however this performs two point locations on: (i) the intial points set and; (ii) the (intial + recieved) point set */
278: PetscScalar *LA_coor;
279: PetscInt bs;
280: DataField PField;
282: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
283: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);
284: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
286: VecDestroy(&pos);
287: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
289: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
290: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
291: for (p=0; p<npoints2; p++) {
292: rankval[p] = LA_sfcell[p].index;
293: }
294: PetscSFDestroy(&sfcell);
295: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
297: /* remove points which left processor */
298: DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
299: DataFieldGetEntries(PField,(void**)&rankval);
301: DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
302: for (p=0; p<npoints2; p++) {
303: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
304: /* kill point */
305: DataBucketRemovePointAtIndex(swarm->db,p);
306: DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
307: DataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
308: p--; /* check replacement point */
309: }
310: }
311: }
312: #endif
314: { /* this performs two point locations: (i) on the intial points set prior to communication; and (ii) on the new (recieved) points */
315: PetscScalar *LA_coor;
316: PetscInt npoints_from_neighbours,bs;
317: DataField PField;
318:
319: npoints_from_neighbours = npoints2 - npoints_prior_migration;
320:
321: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
322: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints_from_neighbours,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);
324: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
325:
326: VecDestroy(&pos);
327: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
328:
329: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
330: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
331: for (p=0; p<npoints_from_neighbours; p++) {
332: rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
333: }
334: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
335: PetscSFDestroy(&sfcell);
336:
337: /* remove points which left processor */
338: DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
339: DataFieldGetEntries(PField,(void**)&rankval);
340:
341: DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
342: for (p=npoints_prior_migration; p<npoints2; p++) {
343: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
344: /* kill point */
345: DataBucketRemovePointAtIndex(swarm->db,p);
346: DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
347: DataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
348: p--; /* check replacement point */
349: }
350: }
351: }
352:
353: {
354: PetscInt *p_cellid;
355:
356: DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
357: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
358: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
359: for (p=0; p<npoints2; p++) {
360: p_cellid[p] = rankval[p];
361: }
362: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
363: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
364: }
365:
366: /* check for error on removed points */
367: if (error_check) {
368: DMSwarmGetSize(dm,&npoints2g);
369: if (npointsg != npoints2g) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points from the DMSwarm must remain constant during migration (initial %D - final %D)",npointsg,npoints2g);
370: }
371: return(0);
372: }
374: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points)
375: {
377: return(0);
378: }
380: /*
381: Redundant as this assumes points can only be sent to a single rank
382: */
383: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
384: {
385: DM_Swarm *swarm = (DM_Swarm*)dm->data;
387: DataEx de;
388: PetscInt p,npoints,*rankval,n_points_recv;
389: PetscMPIInt rank,nrank,negrank;
390: void *point_buffer,*recv_points;
391: size_t sizeof_dmswarm_point;
394: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
395: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
396: *globalsize = npoints;
397: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
398: DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
399: DataExTopologyInitialize(de);
400: for (p=0; p<npoints; p++) {
401: negrank = rankval[p];
402: if (negrank < 0) {
403: nrank = -negrank - 1;
404: DataExTopologyAddNeighbour(de,nrank);
405: }
406: }
407: DataExTopologyFinalize(de);
408: DataExInitializeSendCount(de);
409: for (p=0; p<npoints; p++) {
410: negrank = rankval[p];
411: if (negrank < 0) {
412: nrank = -negrank - 1;
413: DataExAddToSendCount(de,nrank,1);
414: }
415: }
416: DataExFinalizeSendCount(de);
417: DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
418: DataExPackInitialize(de,sizeof_dmswarm_point);
419: for (p=0; p<npoints; p++) {
420: negrank = rankval[p];
421: if (negrank < 0) {
422: nrank = -negrank - 1;
423: rankval[p] = nrank;
424: /* copy point into buffer */
425: DataBucketFillPackedArray(swarm->db,p,point_buffer);
426: /* insert point buffer into DataExchanger */
427: DataExPackData(de,nrank,1,point_buffer);
428: rankval[p] = negrank;
429: }
430: }
431: DataExPackFinalize(de);
432: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
433: DataExBegin(de);
434: DataExEnd(de);
435: DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
436: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
437: DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);
438: for (p=0; p<n_points_recv; p++) {
439: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
441: DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
442: }
443: DataExView(de);
444: DataBucketDestroyPackedArray(swarm->db,&point_buffer);
445: DataExDestroy(de);
446: return(0);
447: }
449: typedef struct {
450: PetscMPIInt owner_rank;
451: PetscReal min[3],max[3];
452: } CollectBBox;
454: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
455: {
456: DM_Swarm *swarm = (DM_Swarm*)dm->data;
458: DataEx de;
459: PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
460: PetscMPIInt rank,nrank;
461: void *point_buffer,*recv_points;
462: size_t sizeof_dmswarm_point,sizeof_bbox_ctx;
463: PetscBool isdmda;
464: CollectBBox *bbox,*recv_bbox;
465: const PetscMPIInt *dmneighborranks;
466: DM dmcell;
469: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
471: DMSwarmGetCellDM(dm,&dmcell);
472: if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
473: isdmda = PETSC_FALSE;
474: PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);
475: if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox");
477: DMGetDimension(dm,&dim);
478: sizeof_bbox_ctx = sizeof(CollectBBox);
479: PetscMalloc1(1,&bbox);
480: bbox->owner_rank = rank;
482: /* compute the bounding box based on the overlapping / stenctil size */
483: {
484: Vec lcoor;
486: DMGetCoordinatesLocal(dmcell,&lcoor);
487: if (dim >= 1) {
488: VecStrideMin(lcoor,0,NULL,&bbox->min[0]);
489: VecStrideMax(lcoor,0,NULL,&bbox->max[0]);
490: }
491: if (dim >= 2) {
492: VecStrideMin(lcoor,1,NULL,&bbox->min[1]);
493: VecStrideMax(lcoor,1,NULL,&bbox->max[1]);
494: }
495: if (dim == 3) {
496: VecStrideMin(lcoor,2,NULL,&bbox->min[2]);
497: VecStrideMax(lcoor,2,NULL,&bbox->max[2]);
498: }
499: }
500: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
501: *globalsize = npoints;
502: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
503: DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
504: /* use DMDA neighbours */
505: DMDAGetNeighbors(dmcell,&dmneighborranks);
506: if (dim == 1) {
507: neighbour_cells = 3;
508: } else if (dim == 2) {
509: neighbour_cells = 9;
510: } else {
511: neighbour_cells = 27;
512: }
513: DataExTopologyInitialize(de);
514: for (p=0; p<neighbour_cells; p++) {
515: if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
516: DataExTopologyAddNeighbour(de,dmneighborranks[p]);
517: }
518: }
519: DataExTopologyFinalize(de);
520: DataExInitializeSendCount(de);
521: for (p=0; p<neighbour_cells; p++) {
522: if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
523: DataExAddToSendCount(de,dmneighborranks[p],1);
524: }
525: }
526: DataExFinalizeSendCount(de);
527: /* send bounding boxes */
528: DataExPackInitialize(de,sizeof_bbox_ctx);
529: for (p=0; p<neighbour_cells; p++) {
530: nrank = dmneighborranks[p];
531: if ( (nrank >= 0) && (nrank != rank) ) {
532: /* insert bbox buffer into DataExchanger */
533: DataExPackData(de,nrank,1,bbox);
534: }
535: }
536: DataExPackFinalize(de);
537: /* recv bounding boxes */
538: DataExBegin(de);
539: DataExEnd(de);
540: DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);
541: for (p=0; p<n_bbox_recv; p++) {
542: PetscPrintf(PETSC_COMM_SELF,"[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank,
543: (double)recv_bbox[p].min[0],(double)recv_bbox[p].max[0],(double)recv_bbox[p].min[1],(double)recv_bbox[p].max[1]);
544: }
545: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
546: DataExInitializeSendCount(de);
547: for (pk=0; pk<n_bbox_recv; pk++) {
548: PetscReal *array_x,*array_y;
550: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
551: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
552: for (p=0; p<npoints; p++) {
553: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
554: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
555: DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);
556: }
557: }
558: }
559: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
560: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
561: }
562: DataExFinalizeSendCount(de);
563: DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
564: DataExPackInitialize(de,sizeof_dmswarm_point);
565: for (pk=0; pk<n_bbox_recv; pk++) {
566: PetscReal *array_x,*array_y;
568: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
569: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
570: for (p=0; p<npoints; p++) {
571: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
572: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
573: /* copy point into buffer */
574: DataBucketFillPackedArray(swarm->db,p,point_buffer);
575: /* insert point buffer into DataExchanger */
576: DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);
577: }
578: }
579: }
580: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
581: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
582: }
583: DataExPackFinalize(de);
584: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
585: DataExBegin(de);
586: DataExEnd(de);
587: DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
588: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
589: DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);
590: for (p=0; p<n_points_recv; p++) {
591: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
593: DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
594: }
595: DataBucketDestroyPackedArray(swarm->db,&point_buffer);
596: PetscFree(bbox);
597: DataExView(de);
598: DataExDestroy(de);
599: return(0);
600: }
603: /* General collection when no order, or neighbour information is provided */
604: /*
605: User provides context and collect() method
606: Broadcast user context
608: for each context / rank {
609: collect(swarm,context,n,list)
610: }
611: */
612: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
613: {
614: DM_Swarm *swarm = (DM_Swarm*)dm->data;
616: DataEx de;
617: PetscInt p,r,npoints,n_points_recv;
618: PetscMPIInt commsize,rank;
619: void *point_buffer,*recv_points;
620: void *ctxlist;
621: PetscInt *n2collect,**collectlist;
622: size_t sizeof_dmswarm_point;
625: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);
626: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
627: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
628: *globalsize = npoints;
629: /* Broadcast user context */
630: PetscMalloc(ctx_size*commsize,&ctxlist);
631: MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));
632: PetscMalloc1(commsize,&n2collect);
633: PetscMalloc1(commsize,&collectlist);
634: for (r=0; r<commsize; r++) {
635: PetscInt _n2collect;
636: PetscInt *_collectlist;
637: void *_ctx_r;
639: _n2collect = 0;
640: _collectlist = NULL;
641: if (r != rank) { /* don't collect data from yourself */
642: _ctx_r = (void*)( (char*)ctxlist + r * ctx_size );
643: collect(dm,_ctx_r,&_n2collect,&_collectlist);
644: }
645: n2collect[r] = _n2collect;
646: collectlist[r] = _collectlist;
647: }
648: DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
649: /* Define topology */
650: DataExTopologyInitialize(de);
651: for (r=0; r<commsize; r++) {
652: if (n2collect[r] > 0) {
653: DataExTopologyAddNeighbour(de,(PetscMPIInt)r);
654: }
655: }
656: DataExTopologyFinalize(de);
657: /* Define send counts */
658: DataExInitializeSendCount(de);
659: for (r=0; r<commsize; r++) {
660: if (n2collect[r] > 0) {
661: DataExAddToSendCount(de,r,n2collect[r]);
662: }
663: }
664: DataExFinalizeSendCount(de);
665: /* Pack data */
666: DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
667: DataExPackInitialize(de,sizeof_dmswarm_point);
668: for (r=0; r<commsize; r++) {
669: for (p=0; p<n2collect[r]; p++) {
670: DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);
671: /* insert point buffer into the data exchanger */
672: DataExPackData(de,r,1,point_buffer);
673: }
674: }
675: DataExPackFinalize(de);
676: /* Scatter */
677: DataExBegin(de);
678: DataExEnd(de);
679: /* Collect data in DMSwarm container */
680: DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
681: DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
682: DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);
683: for (p=0; p<n_points_recv; p++) {
684: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
686: DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
687: }
688: /* Release memory */
689: for (r=0; r<commsize; r++) {
690: if (collectlist[r]) PetscFree(collectlist[r]);
691: }
692: PetscFree(collectlist);
693: PetscFree(n2collect);
694: PetscFree(ctxlist);
695: DataBucketDestroyPackedArray(swarm->db,&point_buffer);
696: DataExView(de);
697: DataExDestroy(de);
698: return(0);
699: }