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