Actual source code: swarmpic_sort.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petscdm.h>
  2:  #include <petscdmda.h>
  3:  #include <petscdmplex.h>
  4:  #include <petscdmswarm.h>
  5:  #include <petsc/private/dmswarmimpl.h>

  7: int sort_CompareSwarmPoint(const void *dataA,const void *dataB)
  8: {
  9:   SwarmPoint *pointA;
 10:   SwarmPoint *pointB;
 11: 
 12:   pointA = (SwarmPoint*)dataA;
 13:   pointB = (SwarmPoint*)dataB;
 14: 
 15:   if (pointA->cell_index < pointB->cell_index) {
 16:     return -1;
 17:   } else if (pointA->cell_index > pointB->cell_index) {
 18:     return 1;
 19:   } else {
 20:     return 0;
 21:   }
 22: }

 24: PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)
 25: {
 27:   qsort(ctx->list,ctx->npoints,sizeof(SwarmPoint),sort_CompareSwarmPoint);
 28:   return(0);
 29: }

 31: PetscErrorCode DMSwarmSortCreate(DMSwarmSort *_ctx)
 32: {
 34:   DMSwarmSort    ctx;
 35: 
 37:   PetscMalloc1(1,&ctx);
 38:   PetscMemzero(ctx,sizeof(struct _p_DMSwarmSort));
 39:   ctx->isvalid = PETSC_FALSE;
 40:   ctx->ncells = 0;
 41:   ctx->npoints = 0;
 42:   PetscMalloc1(1,&ctx->pcell_offsets);
 43:   PetscMalloc1(1,&ctx->list);
 44:   *_ctx = ctx;
 45:   return(0);
 46: }

 48: PetscErrorCode DMSwarmSortSetup(DMSwarmSort ctx,DM dm,PetscInt ncells)
 49: {
 50:   PetscInt        *swarm_cellid;
 51:   PetscInt        p,npoints;
 52:   PetscInt        tmp,c,count;
 53:   PetscErrorCode  ierr;
 54: 
 56:   if (!ctx) return(0);
 57:   if (ctx->isvalid) return(0);
 58: 
 59:   PetscLogEventBegin(DMSWARM_Sort,0,0,0,0);
 60:   /* check the number of cells */
 61:   if (ncells != ctx->ncells) {
 62:     PetscRealloc(sizeof(PetscInt)*(ncells + 1),&ctx->pcell_offsets);
 63:     ctx->ncells = ncells;
 64:   }
 65:   PetscMemzero(ctx->pcell_offsets,sizeof(PetscInt)*(ctx->ncells + 1));
 66: 
 67:   /* get the number of points */
 68:   DMSwarmGetLocalSize(dm,&npoints);
 69:   if (npoints != ctx->npoints) {
 70:     PetscRealloc(sizeof(SwarmPoint)*npoints,&ctx->list);
 71:     ctx->npoints = npoints;
 72:   }
 73:   PetscMemzero(ctx->list,sizeof(SwarmPoint)*npoints);
 74: 
 75:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
 76:   for (p=0; p<ctx->npoints; p++) {
 77:     ctx->list[p].point_index = p;
 78:     ctx->list[p].cell_index  = swarm_cellid[p];
 79:   }
 80:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
 81: 
 82:   DMSwarmSortApplyCellIndexSort(ctx);
 83: 
 84:   /* sum points per cell */
 85:   for (p=0; p<ctx->npoints; p++) {
 86:     ctx->pcell_offsets[ ctx->list[p].cell_index ]++;
 87:   }
 88: 
 89:   /* create offset list */
 90:   count = 0;
 91:   for (c=0; c<ctx->ncells; c++) {
 92:     tmp = ctx->pcell_offsets[c];
 93:     ctx->pcell_offsets[c] = count;
 94:     count = count + tmp;
 95:   }
 96:   ctx->pcell_offsets[c] = count;
 97: 
 98:   ctx->isvalid = PETSC_TRUE;
 99:   PetscLogEventEnd(DMSWARM_Sort,0,0,0,0);
100:   return(0);
101: }

103: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
104: {
105:   DMSwarmSort     ctx;
106:   PetscErrorCode  ierr;
107: 
109:   if (!_ctx) return(0);
110:   if (!*_ctx) return(0);
111:   ctx = *_ctx;
112:   if (ctx->list)      {
113:     PetscFree(ctx->list);
114:   }
115:   if (ctx->pcell_offsets) {
116:     PetscFree(ctx->pcell_offsets);
117:   }
118:   PetscFree(ctx);
119:   *_ctx = NULL;
120:   return(0);
121: }

123: /*@C
124:    DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
125:  
126:    Not collective
127:  
128:    Input parameters:
129: +  dm - a DMSwarm objects
130: .  e - the index of the cell
131: -  npoints - the number of points in the cell
132:  
133:    Level: advanced
134:  
135:    Notes:
136:    You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetNumberOfPointsPerCell()
137:  
138: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess(), DMSwarmSortGetPointsPerCell()
139: @*/
140: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM dm,PetscInt e,PetscInt *npoints)
141: {
142:   DM_Swarm     *swarm = (DM_Swarm*)dm->data;
143:   PetscInt     points_per_cell;
144:   DMSwarmSort  ctx;
145: 
147:   ctx = swarm->sort_context;
148:   if (!ctx) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
149:   if (!ctx->isvalid) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first");
150:   if (e >= ctx->ncells) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cell index (%D) is greater than max number of local cells (%D)",e,ctx->ncells);
151:   if (e < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cell index (%D) cannot be negative",e);
152:   points_per_cell = ctx->pcell_offsets[e+1] - ctx->pcell_offsets[e];
153:   *npoints = points_per_cell;
154:   return(0);
155: }

157: /*@C
158:    DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
159:  
160:    Not collective
161:  
162:    Input parameters:
163: +  dm - a DMSwarm object
164: .  e - the index of the cell
165: .  npoints - the number of points in the cell
166: -  pidlist - array of the indices indentifying all points in cell e
167:  
168:    Level: advanced
169:  
170:    Notes:
171:    - You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell()
172:    - The array pidlist is internally created and must be free'd by the user
173:  
174: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess(), DMSwarmSortGetNumberOfPointsPerCell()
175: @*/
176: PETSC_EXTERN PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm,PetscInt e,PetscInt *npoints,PetscInt **pidlist)
177: {
178:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
180:   PetscInt       points_per_cell;
181:   PetscInt       p,pid,pid_unsorted;
182:   PetscInt       *plist;
183:   DMSwarmSort    ctx;
184: 
186:   ctx = swarm->sort_context;
187:   if (!ctx) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
188:   DMSwarmSortGetNumberOfPointsPerCell(dm,e,&points_per_cell);
189:   PetscMalloc1(points_per_cell,&plist);
190:   for (p=0; p<points_per_cell; p++) {
191:     pid = ctx->pcell_offsets[e] + p;
192:     pid_unsorted = ctx->list[pid].point_index;
193:     plist[p] = pid_unsorted;
194:   }
195:   *npoints = points_per_cell;
196:   *pidlist = plist;
197: 
198:   return(0);
199: }

201: /*@C
202:    DMSwarmSortGetAccess - Setups up a DMSwarm point sort context for efficient traversal of points within a cell
203:  
204:    Not collective
205:  
206:    Input parameter:
207: .  dm - a DMSwarm object
208:  
209:    Calling DMSwarmSortGetAccess() creates a list which enables easy identification of all points contained in a 
210:    given cell. This method does not explicitly sort the data within the DMSwarm based on the cell index associated
211:    with a DMSwarm point. 
212:  
213:    The sort context is valid only for the DMSwarm points defined at the time when DMSwarmSortGetAccess() was called.
214:    For example, suppose the swarm contained NP points when DMSwarmSortGetAccess() was called. If the user subsequently
215:    adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points. 
216:    The indices associated with the 10 new additional points will not be contained within the sort context. 
217:    This means that the user can still safely perform queries via DMSwarmSortGetPointsPerCell() and
218:    DMSwarmSortGetPointsPerCell(), however the results return will be based on the first NP points.
219:  
220:    If any DMSwam re-sizing method is called after DMSwarmSortGetAccess() which modifies any of the first NP entries 
221:    in the DMSwarm, the sort context will become invalid. Currently there are no guards to prevent the user from 
222:    invalidating the sort context. For this reason, we highly recommend you do not use DMSwarmRemovePointAtIndex() in
223:    between calls to DMSwarmSortGetAccess() and DMSwarmSortRestoreAccess().
224:  
225:    To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the 
226:    first pass "marks" points for removal, and the second pass actually removes the points from the DMSwarm.
227:  
228:    Notes:
229:    - You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell() or DMSwarmSortGetNumberOfPointsPerCell()
230:    - The sort context may become invalid if any re-sizing methods are applied which alter the first NP points 
231:      within swarm at the time DMSwarmSortGetAccess() was called.
232:    - You must call DMSwarmSortRestoreAccess() when you no longer need access to the sort context

234:    Level: advanced
235:  
236: .seealso: DMSwarmSetType(), DMSwarmSortRestoreAccess()
237: @*/
238: PETSC_EXTERN PetscErrorCode DMSwarmSortGetAccess(DM dm)
239: {
240:   DM_Swarm        *swarm = (DM_Swarm*)dm->data;
241:   PetscErrorCode  ierr;
242:   PetscInt        ncells;
243:   DM              celldm;
244:   PetscBool       isda,isplex,isshell;
245: 
247:   if (!swarm->sort_context) {
248:     DMSwarmSortCreate(&swarm->sort_context);
249:   }
250: 
251:   /* get the number of cells */
252:   DMSwarmGetCellDM(dm,&celldm);
253:   PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);
254:   PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);
255:   PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);
256:   ncells = 0;
257:   if (isda) {
258:     PetscInt nel,npe;
259:     const PetscInt *element;
260: 
261:     DMDAGetElements(celldm,&nel,&npe,&element);
262:     ncells = nel;
263:     DMDARestoreElements(celldm,&nel,&npe,&element);
264:   } else if (isplex) {
265:     PetscInt ps,pe;
266: 
267:     DMPlexGetHeightStratum(celldm,0,&ps,&pe);
268:     ncells = pe - ps;
269:   } else if (isshell) {
270:     PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);
271: 
272:     PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);
273:     if (method_DMShellGetNumberOfCells) {
274:       method_DMShellGetNumberOfCells(celldm,&ncells);
275:     } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for the DMSHELL object. User must provide a method via PetscObjectComposeFunction( (PetscObject)shelldm, \"DMGetNumberOfCells_C\", your_function_to_compute_number_of_cells );");
276:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
277: 
278:   /* setup */
279:   DMSwarmSortSetup(swarm->sort_context,dm,ncells);
280: 
281:   return(0);
282: }

284: /*@C
285:    DMSwarmSortRestoreAccess - Invalidates the DMSwarm point sorting context
286:  
287:    Not collective
288:  
289:    Input parameter:
290: .  dm - a DMSwarm object
291:  
292:    Level: advanced
293:  
294:    Note:
295:    You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()
296:  
297: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
298: @*/
299: PETSC_EXTERN PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
300: {
301:   DM_Swarm *swarm = (DM_Swarm*)dm->data;
302: 
304:   if (!swarm->sort_context) return(0);
305:   if (!swarm->sort_context->isvalid) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()");
306:   swarm->sort_context->isvalid = PETSC_FALSE;
307:   return(0);
308: }

310: /*@C
311:    DMSwarmSortGetIsValid - Gets the isvalid flag associated with a DMSwarm point sorting context
312:  
313:    Not collective
314:  
315:    Input parameter:
316: .  dm - a DMSwarm object

318:    Output parameter:
319: .  isvalid - flag indicating whether the sort context is up-to-date

321:  Level: advanced
322:  
323: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
324: @*/
325: PETSC_EXTERN PetscErrorCode DMSwarmSortGetIsValid(DM dm,PetscBool *isvalid)
326: {
327:   DM_Swarm *swarm = (DM_Swarm*)dm->data;
328: 
330:   if (!swarm->sort_context) {
331:     *isvalid = PETSC_FALSE;
332:     return(0);
333:   }
334:   *isvalid = swarm->sort_context->isvalid;
335:   return(0);
336: }

338: /*@C
339:    DMSwarmSortGetSizes - Gets the sizes associated with a DMSwarm point sorting context
340:  
341:    Not collective
342:  
343:    Input parameter:
344: .  dm - a DMSwarm object
345:  
346:    Output parameters:
347: +  ncells - number of cells within the sort context (pass NULL to ignore)
348: -  npoints - number of points used to create the sort context (pass NULL to ignore)
349:  
350:    Level: advanced
351:  
352: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
353: @*/
354: PETSC_EXTERN PetscErrorCode DMSwarmSortGetSizes(DM dm,PetscInt *ncells,PetscInt *npoints)
355: {
356:   DM_Swarm *swarm = (DM_Swarm*)dm->data;
357: 
359:   if (!swarm->sort_context) {
360:     if (ncells) { *ncells = 0; }
361:     if (npoints) { *npoints = 0; }
362:     return(0);
363:   }
364:   if (ncells) { *ncells = swarm->sort_context->ncells; }
365:   if (npoints) { *npoints = swarm->sort_context->npoints; }
366:   return(0);
367: }