Actual source code: swarm_migrate.c
petsc-3.14.6 2021-03-30
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;
211: DMSwarmDataBucketGetSizes(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
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: DMSwarmDataBucketGetSizes(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 (size > 1) {
250: DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);
251: } else {
252: DMSwarmDataField PField;
253: PetscInt npoints_curr;
255: /* remove points which the domain */
256: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
257: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
259: DMSwarmDataBucketGetSizes(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: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
264: DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL); /* you need to update npoints as the list size decreases! */
265: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
266: p--; /* check replacement point */
267: }
268: }
269: DMSwarmGetSize(dm,&npoints_prior_migration);
271: }
273: /* locate points newly received */
274: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
276: #if 0
277: { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
278: PetscScalar *LA_coor;
279: PetscInt bs;
280: DMSwarmDataField 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: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
299: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
301: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
302: for (p=0; p<npoints2; p++) {
303: if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
304: /* kill point */
305: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
306: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
307: DMSwarmDataFieldGetEntries(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 initial points set prior to communication; and (ii) on the new (received) points */
315: PetscScalar *LA_coor;
316: PetscInt npoints_from_neighbours,bs;
317: DMSwarmDataField PField;
319: npoints_from_neighbours = npoints2 - npoints_prior_migration;
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);
326: VecDestroy(&pos);
327: DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
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);
337: /* remove points which left processor */
338: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
339: DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
341: DMSwarmDataBucketGetSizes(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: DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
346: DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
347: DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
348: p--; /* check replacement point */
349: }
350: }
351: }
353: {
354: PetscInt *p_cellid;
356: DMSwarmDataBucketGetSizes(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: }
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: DMSwarmDataEx 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: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
396: *globalsize = npoints;
397: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
398: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
399: DMSwarmDataExTopologyInitialize(de);
400: for (p=0; p<npoints; p++) {
401: negrank = rankval[p];
402: if (negrank < 0) {
403: nrank = -negrank - 1;
404: DMSwarmDataExTopologyAddNeighbour(de,nrank);
405: }
406: }
407: DMSwarmDataExTopologyFinalize(de);
408: DMSwarmDataExInitializeSendCount(de);
409: for (p=0; p<npoints; p++) {
410: negrank = rankval[p];
411: if (negrank < 0) {
412: nrank = -negrank - 1;
413: DMSwarmDataExAddToSendCount(de,nrank,1);
414: }
415: }
416: DMSwarmDataExFinalizeSendCount(de);
417: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
418: DMSwarmDataExPackInitialize(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: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
426: /* insert point buffer into DMSwarmDataExchanger */
427: DMSwarmDataExPackData(de,nrank,1,point_buffer);
428: rankval[p] = negrank;
429: }
430: }
431: DMSwarmDataExPackFinalize(de);
432: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
433: DMSwarmDataExBegin(de);
434: DMSwarmDataExEnd(de);
435: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
436: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
437: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_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: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
442: }
443: DMSwarmDataExView(de);
444: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
445: DMSwarmDataExDestroy(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;
457: PetscErrorCode ierr;
458: DMSwarmDataEx 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: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
501: *globalsize = npoints;
502: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
503: DMSwarmDataExCreate(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: DMSwarmDataExTopologyInitialize(de);
514: for (p=0; p<neighbour_cells; p++) {
515: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) {
516: DMSwarmDataExTopologyAddNeighbour(de,dmneighborranks[p]);
517: }
518: }
519: DMSwarmDataExTopologyFinalize(de);
520: DMSwarmDataExInitializeSendCount(de);
521: for (p=0; p<neighbour_cells; p++) {
522: if ((dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank)) {
523: DMSwarmDataExAddToSendCount(de,dmneighborranks[p],1);
524: }
525: }
526: DMSwarmDataExFinalizeSendCount(de);
527: /* send bounding boxes */
528: DMSwarmDataExPackInitialize(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 DMSwarmDataExchanger */
533: DMSwarmDataExPackData(de,nrank,1,bbox);
534: }
535: }
536: DMSwarmDataExPackFinalize(de);
537: /* recv bounding boxes */
538: DMSwarmDataExBegin(de);
539: DMSwarmDataExEnd(de);
540: DMSwarmDataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);
541: /* Wrong, should not be using PETSC_COMM_WORLD */
542: for (p=0; p<n_bbox_recv; p++) {
543: 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,
544: (double)recv_bbox[p].min[0],(double)recv_bbox[p].max[0],(double)recv_bbox[p].min[1],(double)recv_bbox[p].max[1]);
545: }
546: PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
547: /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
548: DMSwarmDataExInitializeSendCount(de);
549: for (pk=0; pk<n_bbox_recv; pk++) {
550: PetscReal *array_x,*array_y;
552: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
553: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
554: for (p=0; p<npoints; p++) {
555: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
556: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
557: DMSwarmDataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);
558: }
559: }
560: }
561: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
562: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
563: }
564: DMSwarmDataExFinalizeSendCount(de);
565: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
566: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
567: for (pk=0; pk<n_bbox_recv; pk++) {
568: PetscReal *array_x,*array_y;
570: DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
571: DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
572: for (p=0; p<npoints; p++) {
573: if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
574: if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
575: /* copy point into buffer */
576: DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
577: /* insert point buffer into DMSwarmDataExchanger */
578: DMSwarmDataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);
579: }
580: }
581: }
582: DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
583: DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
584: }
585: DMSwarmDataExPackFinalize(de);
586: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
587: DMSwarmDataExBegin(de);
588: DMSwarmDataExEnd(de);
589: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
590: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
591: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
592: for (p=0; p<n_points_recv; p++) {
593: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
595: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
596: }
597: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
598: PetscFree(bbox);
599: DMSwarmDataExView(de);
600: DMSwarmDataExDestroy(de);
601: return(0);
602: }
605: /* General collection when no order, or neighbour information is provided */
606: /*
607: User provides context and collect() method
608: Broadcast user context
610: for each context / rank {
611: collect(swarm,context,n,list)
612: }
613: */
614: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
615: {
616: DM_Swarm *swarm = (DM_Swarm*)dm->data;
618: DMSwarmDataEx de;
619: PetscInt p,r,npoints,n_points_recv;
620: PetscMPIInt size,rank;
621: void *point_buffer,*recv_points;
622: void *ctxlist;
623: PetscInt *n2collect,**collectlist;
624: size_t sizeof_dmswarm_point;
627: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
628: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
629: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
630: *globalsize = npoints;
631: /* Broadcast user context */
632: PetscMalloc(ctx_size*size,&ctxlist);
633: MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));
634: PetscMalloc1(size,&n2collect);
635: PetscMalloc1(size,&collectlist);
636: for (r=0; r<size; r++) {
637: PetscInt _n2collect;
638: PetscInt *_collectlist;
639: void *_ctx_r;
641: _n2collect = 0;
642: _collectlist = NULL;
643: if (r != rank) { /* don't collect data from yourself */
644: _ctx_r = (void*)( (char*)ctxlist + r * ctx_size);
645: collect(dm,_ctx_r,&_n2collect,&_collectlist);
646: }
647: n2collect[r] = _n2collect;
648: collectlist[r] = _collectlist;
649: }
650: DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
651: /* Define topology */
652: DMSwarmDataExTopologyInitialize(de);
653: for (r=0; r<size; r++) {
654: if (n2collect[r] > 0) {
655: DMSwarmDataExTopologyAddNeighbour(de,(PetscMPIInt)r);
656: }
657: }
658: DMSwarmDataExTopologyFinalize(de);
659: /* Define send counts */
660: DMSwarmDataExInitializeSendCount(de);
661: for (r=0; r<size; r++) {
662: if (n2collect[r] > 0) {
663: DMSwarmDataExAddToSendCount(de,r,n2collect[r]);
664: }
665: }
666: DMSwarmDataExFinalizeSendCount(de);
667: /* Pack data */
668: DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
669: DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
670: for (r=0; r<size; r++) {
671: for (p=0; p<n2collect[r]; p++) {
672: DMSwarmDataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);
673: /* insert point buffer into the data exchanger */
674: DMSwarmDataExPackData(de,r,1,point_buffer);
675: }
676: }
677: DMSwarmDataExPackFinalize(de);
678: /* Scatter */
679: DMSwarmDataExBegin(de);
680: DMSwarmDataExEnd(de);
681: /* Collect data in DMSwarm container */
682: DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
683: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
684: DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
685: for (p=0; p<n_points_recv; p++) {
686: void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);
688: DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
689: }
690: /* Release memory */
691: for (r=0; r<size; r++) {
692: if (collectlist[r]) PetscFree(collectlist[r]);
693: }
694: PetscFree(collectlist);
695: PetscFree(n2collect);
696: PetscFree(ctxlist);
697: DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
698: DMSwarmDataExView(de);
699: DMSwarmDataExDestroy(de);
700: return(0);
701: }