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