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