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;
15: DMSwarmDataEx 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: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
25: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
26: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0, &de);
27: DMSwarmDataExTopologyInitialize(de);
28: for (p = 0; p < npoints; ++p) {
29: nrank = rankval[p];
30: if (nrank != rank) {
31: DMSwarmDataExTopologyAddNeighbour(de,nrank);
32: }
33: }
34: DMSwarmDataExTopologyFinalize(de);
35: DMSwarmDataExInitializeSendCount(de);
36: for (p=0; p<npoints; p++) {
37: nrank = rankval[p];
38: if (nrank != rank) {
39: DMSwarmDataExAddToSendCount(de,nrank,1);
40: }
41: }
42: DMSwarmDataExFinalizeSendCount(de);
43: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
44: DMSwarmDataExPackInitialize(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: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
50: /* insert point buffer into DMSwarmDataExchanger */
51: DMSwarmDataExPackData(de,nrank,1,point_buffer);
52: }
53: }
54: DMSwarmDataExPackFinalize(de);
55: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
57: if (remove_sent_points) {
58: DMSwarmDataField gfield;
60: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&gfield);
61: DMSwarmDataFieldGetAccess(gfield);
62: DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);
64: /* remove points which left processor */
65: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
66: for (p=0; p<npoints; p++) {
67: nrank = rankval[p];
68: if (nrank != rank) {
69: /* kill point */
70: DMSwarmDataFieldRestoreAccess(gfield);
72: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
74: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
75: DMSwarmDataFieldGetAccess(gfield);
76: DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);
77: p--; /* check replacement point */
78: }
79: }
80: DMSwarmDataFieldRestoreEntries(gfield,(void**)&rankval);
81: DMSwarmDataFieldRestoreAccess(gfield);
82: }
83: DMSwarmDataExBegin(de);
84: DMSwarmDataExEnd(de);
85: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
86: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
87: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_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: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
92: }
93: DMSwarmDataExView(de);
94: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
95: DMSwarmDataExDestroy(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;
102: PetscErrorCode ierr;
103: DMSwarmDataEx 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: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
115: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
116: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
117: DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);
118: DMSwarmDataExTopologyInitialize(de);
119: for (r=0; r<nneighbors; r++) {
120: _rank = neighbourranks[r];
121: if ((_rank != rank) && (_rank > 0)) {
122: DMSwarmDataExTopologyAddNeighbour(de,_rank);
123: }
124: }
125: DMSwarmDataExTopologyFinalize(de);
126: DMSwarmDataExTopologyGetNeighbours(de,&mynneigh,&myneigh);
127: DMSwarmDataExInitializeSendCount(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: DMSwarmDataExAddToSendCount(de,_rank,1);
133: }
134: }
135: }
136: DMSwarmDataExFinalizeSendCount(de);
137: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
138: DMSwarmDataExPackInitialize(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: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
145: /* insert point buffer into DMSwarmDataExchanger */
146: DMSwarmDataExPackData(de,_rank,1,point_buffer);
147: }
148: }
149: }
150: DMSwarmDataExPackFinalize(de);
151: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
152: if (remove_sent_points) {
153: DMSwarmDataField PField;
155: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
156: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
157: /* remove points which left processor */
158: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
159: for (p=0; p<npoints; p++) {
160: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
161: /* kill point */
162: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
163: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
164: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
165: p--; /* check replacement point */
166: }
167: }
168: }
169: DMSwarmDataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);
170: DMSwarmDataExBegin(de);
171: DMSwarmDataExEnd(de);
172: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
173: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
174: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_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: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
179: }
180: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
181: DMSwarmDataExDestroy(de);
182: return(0);
183: }
185: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm,PetscBool remove_sent_points)
186: {
187: DM_Swarm *swarm = (DM_Swarm*)dm->data;
188: PetscErrorCode ierr;
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 size,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),&size);
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: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
211: PetscMalloc1(npoints_curr, &sf_cells);
213: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
214: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
215: for (p=0; p<npoints_curr; p++) {
217: sf_cells[p].rank = 0;
218: sf_cells[p].index = p_cellid[p];
219: if (p_cellid[p] > range) {
220: range = p_cellid[p];
221: }
222: }
223: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
224: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
226: /*PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell);*/
227: PetscSFCreate(PETSC_COMM_SELF,&sfcell);
228: PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
229: }
230: #endif
232: DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
233: DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
234: DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
236: if (error_check) {
237: DMSwarmGetSize(dm,&npointsg);
238: }
239: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
240: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
241: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
242: for (p=0; p<npoints; p++) {
243: rankval[p] = LA_sfcell[p].index;
244: }
245: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
246: PetscSFDestroy(&sfcell);
248: if (size > 1) {
249: DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);
250: } else {
251: DMSwarmDataField PField;
252: PetscInt npoints_curr;
254: /* remove points which the domain */
255: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
256: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
258: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
259: for (p=0; p<npoints_curr; p++) {
260: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
261: /* kill point */
262: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
263: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL); /* you need to update npoints as the list size decreases! */
264: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
265: p--; /* check replacement point */
266: }
267: }
268: DMSwarmGetSize(dm,&npoints_prior_migration);
270: }
272: /* locate points newly received */
273: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
275: #if 0
276: { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
277: PetscScalar *LA_coor;
278: PetscInt bs;
279: DMSwarmDataField PField;
281: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
282: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);
283: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
285: VecDestroy(&pos);
286: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
288: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
289: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
290: for (p=0; p<npoints2; p++) {
291: rankval[p] = LA_sfcell[p].index;
292: }
293: PetscSFDestroy(&sfcell);
294: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
296: /* remove points which left processor */
297: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
298: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
300: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
301: for (p=0; p<npoints2; p++) {
302: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
303: /* kill point */
304: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
305: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
306: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
307: p--; /* check replacement point */
308: }
309: }
310: }
311: #endif
313: { /* this performs two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
314: PetscScalar *LA_coor;
315: PetscInt npoints_from_neighbours,bs;
316: DMSwarmDataField PField;
318: npoints_from_neighbours = npoints2 - npoints_prior_migration;
320: DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
321: VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints_from_neighbours,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);
323: DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
325: VecDestroy(&pos);
326: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
328: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
329: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
330: for (p=0; p<npoints_from_neighbours; p++) {
331: rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
332: }
333: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
334: PetscSFDestroy(&sfcell);
336: /* remove points which left processor */
337: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
338: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
340: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
341: for (p=npoints_prior_migration; p<npoints2; p++) {
342: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
343: /* kill point */
344: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
345: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
346: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
347: p--; /* check replacement point */
348: }
349: }
350: }
352: {
353: PetscInt *p_cellid;
355: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
356: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
357: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
358: for (p=0; p<npoints2; p++) {
359: p_cellid[p] = rankval[p];
360: }
361: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
362: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
363: }
365: /* check for error on removed points */
366: if (error_check) {
367: DMSwarmGetSize(dm,&npoints2g);
368: 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);
369: }
370: return(0);
371: }
373: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points)
374: {
376: return(0);
377: }
379: /*
380: Redundant as this assumes points can only be sent to a single rank
381: */
382: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
383: {
384: DM_Swarm *swarm = (DM_Swarm*)dm->data;
386: DMSwarmDataEx de;
387: PetscInt p,npoints,*rankval,n_points_recv;
388: PetscMPIInt rank,nrank,negrank;
389: void *point_buffer,*recv_points;
390: size_t sizeof_dmswarm_point;
393: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
394: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
395: *globalsize = npoints;
396: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
397: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
398: DMSwarmDataExTopologyInitialize(de);
399: for (p=0; p<npoints; p++) {
400: negrank = rankval[p];
401: if (negrank < 0) {
402: nrank = -negrank - 1;
403: DMSwarmDataExTopologyAddNeighbour(de,nrank);
404: }
405: }
406: DMSwarmDataExTopologyFinalize(de);
407: DMSwarmDataExInitializeSendCount(de);
408: for (p=0; p<npoints; p++) {
409: negrank = rankval[p];
410: if (negrank < 0) {
411: nrank = -negrank - 1;
412: DMSwarmDataExAddToSendCount(de,nrank,1);
413: }
414: }
415: DMSwarmDataExFinalizeSendCount(de);
416: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
417: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
418: for (p=0; p<npoints; p++) {
419: negrank = rankval[p];
420: if (negrank < 0) {
421: nrank = -negrank - 1;
422: rankval[p] = nrank;
423: /* copy point into buffer */
424: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
425: /* insert point buffer into DMSwarmDataExchanger */
426: DMSwarmDataExPackData(de,nrank,1,point_buffer);
427: rankval[p] = negrank;
428: }
429: }
430: DMSwarmDataExPackFinalize(de);
431: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
432: DMSwarmDataExBegin(de);
433: DMSwarmDataExEnd(de);
434: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
435: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
436: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
437: for (p=0; p<n_points_recv; p++) {
438: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
440: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
441: }
442: DMSwarmDataExView(de);
443: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
444: DMSwarmDataExDestroy(de);
445: return(0);
446: }
448: typedef struct {
449: PetscMPIInt owner_rank;
450: PetscReal min[3],max[3];
451: } CollectBBox;
453: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
454: {
455: DM_Swarm * swarm = (DM_Swarm*)dm->data;
456: PetscErrorCode ierr;
457: DMSwarmDataEx de;
458: PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
459: PetscMPIInt rank,nrank;
460: void *point_buffer,*recv_points;
461: size_t sizeof_dmswarm_point,sizeof_bbox_ctx;
462: PetscBool isdmda;
463: CollectBBox *bbox,*recv_bbox;
464: const PetscMPIInt *dmneighborranks;
465: DM dmcell;
468: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
470: DMSwarmGetCellDM(dm,&dmcell);
471: if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
472: isdmda = PETSC_FALSE;
473: PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);
474: if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox");
476: DMGetDimension(dm,&dim);
477: sizeof_bbox_ctx = sizeof(CollectBBox);
478: PetscMalloc1(1,&bbox);
479: bbox->owner_rank = rank;
481: /* compute the bounding box based on the overlapping / stenctil size */
482: {
483: Vec lcoor;
485: DMGetCoordinatesLocal(dmcell,&lcoor);
486: if (dim >= 1) {
487: VecStrideMin(lcoor,0,NULL,&bbox->min[0]);
488: VecStrideMax(lcoor,0,NULL,&bbox->max[0]);
489: }
490: if (dim >= 2) {
491: VecStrideMin(lcoor,1,NULL,&bbox->min[1]);
492: VecStrideMax(lcoor,1,NULL,&bbox->max[1]);
493: }
494: if (dim == 3) {
495: VecStrideMin(lcoor,2,NULL,&bbox->min[2]);
496: VecStrideMax(lcoor,2,NULL,&bbox->max[2]);
497: }
498: }
499: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
500: *globalsize = npoints;
501: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
502: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
503: /* use DMDA neighbours */
504: DMDAGetNeighbors(dmcell,&dmneighborranks);
505: if (dim == 1) {
506: neighbour_cells = 3;
507: } else if (dim == 2) {
508: neighbour_cells = 9;
509: } else {
510: neighbour_cells = 27;
511: }
512: DMSwarmDataExTopologyInitialize(de);
513: for (p=0; p<neighbour_cells; p++) {
514: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) {
515: DMSwarmDataExTopologyAddNeighbour(de,dmneighborranks[p]);
516: }
517: }
518: DMSwarmDataExTopologyFinalize(de);
519: DMSwarmDataExInitializeSendCount(de);
520: for (p=0; p<neighbour_cells; p++) {
521: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) {
522: DMSwarmDataExAddToSendCount(de,dmneighborranks[p],1);
523: }
524: }
525: DMSwarmDataExFinalizeSendCount(de);
526: /* send bounding boxes */
527: DMSwarmDataExPackInitialize(de,sizeof_bbox_ctx);
528: for (p=0; p<neighbour_cells; p++) {
529: nrank = dmneighborranks[p];
530: if ((nrank >= 0) && (nrank != rank)) {
531: /* insert bbox buffer into DMSwarmDataExchanger */
532: DMSwarmDataExPackData(de,nrank,1,bbox);
533: }
534: }
535: DMSwarmDataExPackFinalize(de);
536: /* recv bounding boxes */
537: DMSwarmDataExBegin(de);
538: DMSwarmDataExEnd(de);
539: DMSwarmDataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);
540: /* Wrong, should not be using PETSC_COMM_WORLD */
541: for (p=0; p<n_bbox_recv; p++) {
542: 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,
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: PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
546: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
547: DMSwarmDataExInitializeSendCount(de);
548: for (pk=0; pk<n_bbox_recv; pk++) {
549: PetscReal *array_x,*array_y;
551: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
552: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
553: for (p=0; p<npoints; p++) {
554: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
555: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
556: DMSwarmDataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);
557: }
558: }
559: }
560: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
561: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
562: }
563: DMSwarmDataExFinalizeSendCount(de);
564: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
565: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
566: for (pk=0; pk<n_bbox_recv; pk++) {
567: PetscReal *array_x,*array_y;
569: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
570: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
571: for (p=0; p<npoints; p++) {
572: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
573: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
574: /* copy point into buffer */
575: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
576: /* insert point buffer into DMSwarmDataExchanger */
577: DMSwarmDataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);
578: }
579: }
580: }
581: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
582: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
583: }
584: DMSwarmDataExPackFinalize(de);
585: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
586: DMSwarmDataExBegin(de);
587: DMSwarmDataExEnd(de);
588: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
589: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
590: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
591: for (p=0; p<n_points_recv; p++) {
592: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
594: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
595: }
596: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
597: PetscFree(bbox);
598: DMSwarmDataExView(de);
599: DMSwarmDataExDestroy(de);
600: return(0);
601: }
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: DMSwarmDataEx de;
617: PetscInt p,r,npoints,n_points_recv;
618: PetscMPIInt size,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),&size);
626: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
627: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
628: *globalsize = npoints;
629: /* Broadcast user context */
630: PetscMalloc(ctx_size*size,&ctxlist);
631: MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));
632: PetscMalloc1(size,&n2collect);
633: PetscMalloc1(size,&collectlist);
634: for (r=0; r<size; 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: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
649: /* Define topology */
650: DMSwarmDataExTopologyInitialize(de);
651: for (r=0; r<size; r++) {
652: if (n2collect[r] > 0) {
653: DMSwarmDataExTopologyAddNeighbour(de,(PetscMPIInt)r);
654: }
655: }
656: DMSwarmDataExTopologyFinalize(de);
657: /* Define send counts */
658: DMSwarmDataExInitializeSendCount(de);
659: for (r=0; r<size; r++) {
660: if (n2collect[r] > 0) {
661: DMSwarmDataExAddToSendCount(de,r,n2collect[r]);
662: }
663: }
664: DMSwarmDataExFinalizeSendCount(de);
665: /* Pack data */
666: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
667: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
668: for (r=0; r<size; r++) {
669: for (p=0; p<n2collect[r]; p++) {
670: DMSwarmDataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);
671: /* insert point buffer into the data exchanger */
672: DMSwarmDataExPackData(de,r,1,point_buffer);
673: }
674: }
675: DMSwarmDataExPackFinalize(de);
676: /* Scatter */
677: DMSwarmDataExBegin(de);
678: DMSwarmDataExEnd(de);
679: /* Collect data in DMSwarm container */
680: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
681: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
682: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_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: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
687: }
688: /* Release memory */
689: for (r=0; r<size; r++) {
690: if (collectlist[r]) PetscFree(collectlist[r]);
691: }
692: PetscFree(collectlist);
693: PetscFree(n2collect);
694: PetscFree(ctxlist);
695: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
696: DMSwarmDataExView(de);
697: DMSwarmDataExDestroy(de);
698: return(0);
699: }