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