1 #include <petscdm.h> 2 #include <petscdmda.h> 3 #include <petscdmplex.h> 4 #include <petscdmswarm.h> 5 #include <petsc/private/dmswarmimpl.h> 6 7 static 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 static PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx) 22 { 23 PetscFunctionBegin; 24 qsort(ctx->list, ctx->npoints, sizeof(SwarmPoint), sort_CompareSwarmPoint); 25 PetscFunctionReturn(PETSC_SUCCESS); 26 } 27 28 static 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 static 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: `DMSWARM`, `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: `DMSWARM`, `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 Level: advanced 196 197 Notes: 198 Calling `DMSwarmSortGetAccess()` creates a list which enables easy identification of all points contained in a 199 given cell. This method does not explicitly sort the data within the `DMSWARM` based on the cell index associated 200 with a `DMSWARM` point. 201 202 The sort context is valid only for the `DMSWARM` points defined at the time when `DMSwarmSortGetAccess()` was called. 203 For example, suppose the swarm contained NP points when `DMSwarmSortGetAccess()` was called. If the user subsequently 204 adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points. 205 The indices associated with the 10 new additional points will not be contained within the sort context. 206 This means that the user can still safely perform queries via `DMSwarmSortGetPointsPerCell()` and 207 `DMSwarmSortGetPointsPerCell()`, however the results return will be based on the first NP points. 208 209 If any` DMSWARM` re-sizing method is called after `DMSwarmSortGetAccess()` which modifies any of the first NP entries 210 in the `DMSWARM`, the sort context will become invalid. Currently there are no guards to prevent the user from 211 invalidating the sort context. For this reason, we highly recommend you do not use `DMSwarmRemovePointAtIndex()` in 212 between calls to `DMSwarmSortGetAccess()` and `DMSwarmSortRestoreAccess()`. 213 214 To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the 215 first pass "marks" points for removal, and the second pass actually removes the points from the `DMSWARM` 216 217 You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()` or `DMSwarmSortGetNumberOfPointsPerCell()` 218 219 The sort context may become invalid if any re-sizing methods are applied which alter the first NP points 220 within swarm at the time `DMSwarmSortGetAccess()` was called. 221 222 You must call `DMSwarmSortRestoreAccess()` when you no longer need access to the sort context 223 224 .seealso: `DMSWARM`, `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: `DMSWARM`, `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: `DMSWARM`, `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: `DMSWARM`, `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