1 #include <petscdm.h> 2 #include <petscdmda.h> 3 #include <petscdmplex.h> 4 #include <petscdmswarm.h> 5 #include <petsc/private/dmswarmimpl.h> 6 7 int sort_CompareSwarmPoint(const void *dataA,const void *dataB) 8 { 9 SwarmPoint *pointA = (SwarmPoint*)dataA; 10 SwarmPoint *pointB = (SwarmPoint*)dataB; 11 12 if (pointA->cell_index < pointB->cell_index) { 13 return -1; 14 } else if (pointA->cell_index > pointB->cell_index) { 15 return 1; 16 } else { 17 return 0; 18 } 19 } 20 21 PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx) 22 { 23 PetscFunctionBegin; 24 qsort(ctx->list,ctx->npoints,sizeof(SwarmPoint),sort_CompareSwarmPoint); 25 PetscFunctionReturn(0); 26 } 27 28 PetscErrorCode DMSwarmSortCreate(DMSwarmSort *_ctx) 29 { 30 DMSwarmSort ctx; 31 32 PetscFunctionBegin; 33 PetscCall(PetscNew(&ctx)); 34 ctx->isvalid = PETSC_FALSE; 35 ctx->ncells = 0; 36 ctx->npoints = 0; 37 PetscCall(PetscMalloc1(1,&ctx->pcell_offsets)); 38 PetscCall(PetscMalloc1(1,&ctx->list)); 39 *_ctx = ctx; 40 PetscFunctionReturn(0); 41 } 42 43 PetscErrorCode DMSwarmSortSetup(DMSwarmSort ctx,DM dm,PetscInt ncells) 44 { 45 PetscInt *swarm_cellid; 46 PetscInt p,npoints; 47 PetscInt tmp,c,count; 48 49 PetscFunctionBegin; 50 if (!ctx) PetscFunctionReturn(0); 51 if (ctx->isvalid) PetscFunctionReturn(0); 52 53 PetscCall(PetscLogEventBegin(DMSWARM_Sort,0,0,0,0)); 54 /* check the number of cells */ 55 if (ncells != ctx->ncells) { 56 PetscCall(PetscRealloc(sizeof(PetscInt)*(ncells + 1),&ctx->pcell_offsets)); 57 ctx->ncells = ncells; 58 } 59 PetscCall(PetscArrayzero(ctx->pcell_offsets,ctx->ncells + 1)); 60 61 /* get the number of points */ 62 PetscCall(DMSwarmGetLocalSize(dm,&npoints)); 63 if (npoints != ctx->npoints) { 64 PetscCall(PetscRealloc(sizeof(SwarmPoint)*npoints,&ctx->list)); 65 ctx->npoints = npoints; 66 } 67 PetscCall(PetscArrayzero(ctx->list,npoints)); 68 69 PetscCall(DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid)); 70 for (p=0; p<ctx->npoints; p++) { 71 ctx->list[p].point_index = p; 72 ctx->list[p].cell_index = swarm_cellid[p]; 73 } 74 PetscCall(DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid)); 75 PetscCall(DMSwarmSortApplyCellIndexSort(ctx)); 76 77 /* sum points per cell */ 78 for (p=0; p<ctx->npoints; p++) { 79 ctx->pcell_offsets[ ctx->list[p].cell_index ]++; 80 } 81 82 /* create offset list */ 83 count = 0; 84 for (c=0; c<ctx->ncells; c++) { 85 tmp = ctx->pcell_offsets[c]; 86 ctx->pcell_offsets[c] = count; 87 count = count + tmp; 88 } 89 ctx->pcell_offsets[c] = count; 90 91 ctx->isvalid = PETSC_TRUE; 92 PetscCall(PetscLogEventEnd(DMSWARM_Sort,0,0,0,0)); 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx) 97 { 98 DMSwarmSort ctx; 99 100 PetscFunctionBegin; 101 if (!_ctx) PetscFunctionReturn(0); 102 if (!*_ctx) PetscFunctionReturn(0); 103 ctx = *_ctx; 104 if (ctx->list) { 105 PetscCall(PetscFree(ctx->list)); 106 } 107 if (ctx->pcell_offsets) { 108 PetscCall(PetscFree(ctx->pcell_offsets)); 109 } 110 PetscCall(PetscFree(ctx)); 111 *_ctx = NULL; 112 PetscFunctionReturn(0); 113 } 114 115 /*@C 116 DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell 117 118 Not collective 119 120 Input parameters: 121 + dm - a DMSwarm objects 122 . e - the index of the cell 123 - npoints - the number of points in the cell 124 125 Level: advanced 126 127 Notes: 128 You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetNumberOfPointsPerCell() 129 130 .seealso: DMSwarmSetType(), DMSwarmSortGetAccess(), DMSwarmSortGetPointsPerCell() 131 @*/ 132 PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM dm,PetscInt e,PetscInt *npoints) 133 { 134 DM_Swarm *swarm = (DM_Swarm*)dm->data; 135 PetscInt points_per_cell; 136 DMSwarmSort ctx; 137 138 PetscFunctionBegin; 139 ctx = swarm->sort_context; 140 PetscCheck(ctx,PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first"); 141 PetscCheck(ctx->isvalid,PETSC_COMM_SELF,PETSC_ERR_USER,"SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first"); 142 PetscCheckFalse(e >= ctx->ncells,PETSC_COMM_SELF,PETSC_ERR_USER,"Cell index (%D) is greater than max number of local cells (%D)",e,ctx->ncells); 143 PetscCheckFalse(e < 0,PETSC_COMM_SELF,PETSC_ERR_USER,"Cell index (%D) cannot be negative",e); 144 points_per_cell = ctx->pcell_offsets[e+1] - ctx->pcell_offsets[e]; 145 *npoints = points_per_cell; 146 PetscFunctionReturn(0); 147 } 148 149 /*@C 150 DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell 151 152 Not collective 153 154 Input parameters: 155 + dm - a DMSwarm object 156 . e - the index of the cell 157 . npoints - the number of points in the cell 158 - pidlist - array of the indices indentifying all points in cell e 159 160 Level: advanced 161 162 Notes: 163 You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell() 164 165 The array pidlist is internally created and must be free'd by the user 166 167 .seealso: DMSwarmSetType(), DMSwarmSortGetAccess(), DMSwarmSortGetNumberOfPointsPerCell() 168 @*/ 169 PETSC_EXTERN PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm,PetscInt e,PetscInt *npoints,PetscInt **pidlist) 170 { 171 DM_Swarm *swarm = (DM_Swarm*)dm->data; 172 PetscInt points_per_cell; 173 PetscInt p,pid,pid_unsorted; 174 PetscInt *plist; 175 DMSwarmSort ctx; 176 177 PetscFunctionBegin; 178 ctx = swarm->sort_context; 179 PetscCheck(ctx,PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first"); 180 PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm,e,&points_per_cell)); 181 PetscCall(PetscMalloc1(points_per_cell,&plist)); 182 for (p=0; p<points_per_cell; p++) { 183 pid = ctx->pcell_offsets[e] + p; 184 pid_unsorted = ctx->list[pid].point_index; 185 plist[p] = pid_unsorted; 186 } 187 *npoints = points_per_cell; 188 *pidlist = plist; 189 190 PetscFunctionReturn(0); 191 } 192 193 /*@C 194 DMSwarmSortGetAccess - Setups up a DMSwarm point sort context for efficient traversal of points within a cell 195 196 Not collective 197 198 Input parameter: 199 . dm - a DMSwarm object 200 201 Calling DMSwarmSortGetAccess() creates a list which enables easy identification of all points contained in a 202 given cell. This method does not explicitly sort the data within the DMSwarm based on the cell index associated 203 with a DMSwarm point. 204 205 The sort context is valid only for the DMSwarm points defined at the time when DMSwarmSortGetAccess() was called. 206 For example, suppose the swarm contained NP points when DMSwarmSortGetAccess() was called. If the user subsequently 207 adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points. 208 The indices associated with the 10 new additional points will not be contained within the sort context. 209 This means that the user can still safely perform queries via DMSwarmSortGetPointsPerCell() and 210 DMSwarmSortGetPointsPerCell(), however the results return will be based on the first NP points. 211 212 If any DMSwam re-sizing method is called after DMSwarmSortGetAccess() which modifies any of the first NP entries 213 in the DMSwarm, the sort context will become invalid. Currently there are no guards to prevent the user from 214 invalidating the sort context. For this reason, we highly recommend you do not use DMSwarmRemovePointAtIndex() in 215 between calls to DMSwarmSortGetAccess() and DMSwarmSortRestoreAccess(). 216 217 To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the 218 first pass "marks" points for removal, and the second pass actually removes the points from the DMSwarm. 219 220 Notes: 221 You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell() or DMSwarmSortGetNumberOfPointsPerCell() 222 223 The sort context may become invalid if any re-sizing methods are applied which alter the first NP points 224 within swarm at the time DMSwarmSortGetAccess() was called. 225 226 You must call DMSwarmSortRestoreAccess() when you no longer need access to the sort context 227 228 Level: advanced 229 230 .seealso: DMSwarmSetType(), DMSwarmSortRestoreAccess() 231 @*/ 232 PETSC_EXTERN PetscErrorCode DMSwarmSortGetAccess(DM dm) 233 { 234 DM_Swarm *swarm = (DM_Swarm*)dm->data; 235 PetscInt ncells; 236 DM celldm; 237 PetscBool isda,isplex,isshell; 238 239 PetscFunctionBegin; 240 if (!swarm->sort_context) { 241 PetscCall(DMSwarmSortCreate(&swarm->sort_context)); 242 } 243 244 /* get the number of cells */ 245 PetscCall(DMSwarmGetCellDM(dm,&celldm)); 246 PetscCall(PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda)); 247 PetscCall(PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex)); 248 PetscCall(PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell)); 249 ncells = 0; 250 if (isda) { 251 PetscInt nel,npe; 252 const PetscInt *element; 253 254 PetscCall(DMDAGetElements(celldm,&nel,&npe,&element)); 255 ncells = nel; 256 PetscCall(DMDARestoreElements(celldm,&nel,&npe,&element)); 257 } else if (isplex) { 258 PetscInt ps,pe; 259 260 PetscCall(DMPlexGetHeightStratum(celldm,0,&ps,&pe)); 261 ncells = pe - ps; 262 } else if (isshell) { 263 PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*); 264 265 PetscCall(PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells)); 266 if (method_DMShellGetNumberOfCells) { 267 PetscCall(method_DMShellGetNumberOfCells(celldm,&ncells)); 268 } 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);"); 269 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL"); 270 271 /* setup */ 272 PetscCall(DMSwarmSortSetup(swarm->sort_context,dm,ncells)); 273 PetscFunctionReturn(0); 274 } 275 276 /*@C 277 DMSwarmSortRestoreAccess - Invalidates the DMSwarm point sorting context 278 279 Not collective 280 281 Input parameter: 282 . dm - a DMSwarm object 283 284 Level: advanced 285 286 Note: 287 You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess() 288 289 .seealso: DMSwarmSetType(), DMSwarmSortGetAccess() 290 @*/ 291 PETSC_EXTERN PetscErrorCode DMSwarmSortRestoreAccess(DM dm) 292 { 293 DM_Swarm *swarm = (DM_Swarm*)dm->data; 294 295 PetscFunctionBegin; 296 if (!swarm->sort_context) PetscFunctionReturn(0); 297 PetscCheck(swarm->sort_context->isvalid,PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()"); 298 swarm->sort_context->isvalid = PETSC_FALSE; 299 PetscFunctionReturn(0); 300 } 301 302 /*@C 303 DMSwarmSortGetIsValid - Gets the isvalid flag associated with a DMSwarm point sorting context 304 305 Not collective 306 307 Input parameter: 308 . dm - a DMSwarm object 309 310 Output parameter: 311 . isvalid - flag indicating whether the sort context is up-to-date 312 313 Level: advanced 314 315 .seealso: DMSwarmSetType(), DMSwarmSortGetAccess() 316 @*/ 317 PETSC_EXTERN PetscErrorCode DMSwarmSortGetIsValid(DM dm,PetscBool *isvalid) 318 { 319 DM_Swarm *swarm = (DM_Swarm*)dm->data; 320 321 PetscFunctionBegin; 322 if (!swarm->sort_context) { 323 *isvalid = PETSC_FALSE; 324 PetscFunctionReturn(0); 325 } 326 *isvalid = swarm->sort_context->isvalid; 327 PetscFunctionReturn(0); 328 } 329 330 /*@C 331 DMSwarmSortGetSizes - Gets the sizes associated with a DMSwarm point sorting context 332 333 Not collective 334 335 Input parameter: 336 . dm - a DMSwarm object 337 338 Output parameters: 339 + ncells - number of cells within the sort context (pass NULL to ignore) 340 - npoints - number of points used to create the sort context (pass NULL to ignore) 341 342 Level: advanced 343 344 .seealso: DMSwarmSetType(), DMSwarmSortGetAccess() 345 @*/ 346 PETSC_EXTERN PetscErrorCode DMSwarmSortGetSizes(DM dm,PetscInt *ncells,PetscInt *npoints) 347 { 348 DM_Swarm *swarm = (DM_Swarm*)dm->data; 349 350 PetscFunctionBegin; 351 if (!swarm->sort_context) { 352 if (ncells) *ncells = 0; 353 if (npoints) *npoints = 0; 354 PetscFunctionReturn(0); 355 } 356 if (ncells) *ncells = swarm->sort_context->ncells; 357 if (npoints) *npoints = swarm->sort_context->npoints; 358 PetscFunctionReturn(0); 359 } 360