1 2 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 3 4 static PetscErrorCode DMDAGetElements_1D(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[]) { 5 DM_DA *da = (DM_DA *)dm->data; 6 PetscInt i, xs, xe, Xs, Xe; 7 PetscInt cnt = 0; 8 9 PetscFunctionBegin; 10 if (!da->e) { 11 PetscInt corners[2]; 12 13 PetscCheck(da->s, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot get elements for DMDA with zero stencil width"); 14 PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xe, NULL, NULL)); 15 PetscCall(DMDAGetGhostCorners(dm, &Xs, NULL, NULL, &Xe, NULL, NULL)); 16 xe += xs; 17 Xe += Xs; 18 if (xs != Xs) xs -= 1; 19 da->ne = 1 * (xe - xs - 1); 20 PetscCall(PetscMalloc1(1 + 2 * da->ne, &da->e)); 21 for (i = xs; i < xe - 1; i++) { 22 da->e[cnt++] = (i - Xs); 23 da->e[cnt++] = (i - Xs + 1); 24 } 25 da->nen = 2; 26 27 corners[0] = (xs - Xs); 28 corners[1] = (xe - 1 - Xs); 29 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 2, corners, PETSC_COPY_VALUES, &da->ecorners)); 30 } 31 *nel = da->ne; 32 *nen = da->nen; 33 *e = da->e; 34 PetscFunctionReturn(0); 35 } 36 37 static PetscErrorCode DMDAGetElements_2D(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[]) { 38 DM_DA *da = (DM_DA *)dm->data; 39 PetscInt i, xs, xe, Xs, Xe; 40 PetscInt j, ys, ye, Ys, Ye; 41 PetscInt cnt = 0, cell[4], ns = 2; 42 PetscInt c, split[] = {0, 1, 3, 2, 3, 1}; 43 44 PetscFunctionBegin; 45 if (!da->e) { 46 PetscInt corners[4], nn = 0; 47 48 PetscCheck(da->s, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot get elements for DMDA with zero stencil width"); 49 50 switch (da->elementtype) { 51 case DMDA_ELEMENT_Q1: da->nen = 4; break; 52 case DMDA_ELEMENT_P1: da->nen = 3; break; 53 default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown element type %d", da->elementtype); 54 } 55 nn = da->nen; 56 57 if (da->elementtype == DMDA_ELEMENT_P1) ns = 2; 58 if (da->elementtype == DMDA_ELEMENT_Q1) ns = 1; 59 PetscCall(DMDAGetCorners(dm, &xs, &ys, NULL, &xe, &ye, NULL)); 60 PetscCall(DMDAGetGhostCorners(dm, &Xs, &Ys, NULL, &Xe, &Ye, NULL)); 61 xe += xs; 62 Xe += Xs; 63 if (xs != Xs) xs -= 1; 64 ye += ys; 65 Ye += Ys; 66 if (ys != Ys) ys -= 1; 67 da->ne = ns * (xe - xs - 1) * (ye - ys - 1); 68 PetscCall(PetscMalloc1(1 + nn * da->ne, &da->e)); 69 for (j = ys; j < ye - 1; j++) { 70 for (i = xs; i < xe - 1; i++) { 71 cell[0] = (i - Xs) + (j - Ys) * (Xe - Xs); 72 cell[1] = (i - Xs + 1) + (j - Ys) * (Xe - Xs); 73 cell[2] = (i - Xs + 1) + (j - Ys + 1) * (Xe - Xs); 74 cell[3] = (i - Xs) + (j - Ys + 1) * (Xe - Xs); 75 if (da->elementtype == DMDA_ELEMENT_P1) { 76 for (c = 0; c < ns * nn; c++) da->e[cnt++] = cell[split[c]]; 77 } 78 if (da->elementtype == DMDA_ELEMENT_Q1) { 79 for (c = 0; c < ns * nn; c++) da->e[cnt++] = cell[c]; 80 } 81 } 82 } 83 84 corners[0] = (xs - Xs) + (ys - Ys) * (Xe - Xs); 85 corners[1] = (xe - 1 - Xs) + (ys - Ys) * (Xe - Xs); 86 corners[2] = (xs - Xs) + (ye - 1 - Ys) * (Xe - Xs); 87 corners[3] = (xe - 1 - Xs) + (ye - 1 - Ys) * (Xe - Xs); 88 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 4, corners, PETSC_COPY_VALUES, &da->ecorners)); 89 } 90 *nel = da->ne; 91 *nen = da->nen; 92 *e = da->e; 93 PetscFunctionReturn(0); 94 } 95 96 static PetscErrorCode DMDAGetElements_3D(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[]) { 97 DM_DA *da = (DM_DA *)dm->data; 98 PetscInt i, xs, xe, Xs, Xe; 99 PetscInt j, ys, ye, Ys, Ye; 100 PetscInt k, zs, ze, Zs, Ze; 101 PetscInt cnt = 0, cell[8], ns = 6; 102 PetscInt c, split[] = {0, 1, 3, 7, 0, 1, 7, 4, 1, 2, 3, 7, 1, 2, 7, 6, 1, 4, 5, 7, 1, 5, 6, 7}; 103 104 PetscFunctionBegin; 105 if (!da->e) { 106 PetscInt corners[8], nn = 0; 107 108 PetscCheck(da->s, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot get elements for DMDA with zero stencil width"); 109 110 switch (da->elementtype) { 111 case DMDA_ELEMENT_Q1: da->nen = 8; break; 112 case DMDA_ELEMENT_P1: da->nen = 4; break; 113 default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown element type %d", da->elementtype); 114 } 115 nn = da->nen; 116 117 if (da->elementtype == DMDA_ELEMENT_P1) ns = 6; 118 if (da->elementtype == DMDA_ELEMENT_Q1) ns = 1; 119 PetscCall(DMDAGetCorners(dm, &xs, &ys, &zs, &xe, &ye, &ze)); 120 PetscCall(DMDAGetGhostCorners(dm, &Xs, &Ys, &Zs, &Xe, &Ye, &Ze)); 121 xe += xs; 122 Xe += Xs; 123 if (xs != Xs) xs -= 1; 124 ye += ys; 125 Ye += Ys; 126 if (ys != Ys) ys -= 1; 127 ze += zs; 128 Ze += Zs; 129 if (zs != Zs) zs -= 1; 130 da->ne = ns * (xe - xs - 1) * (ye - ys - 1) * (ze - zs - 1); 131 PetscCall(PetscMalloc1(1 + nn * da->ne, &da->e)); 132 for (k = zs; k < ze - 1; k++) { 133 for (j = ys; j < ye - 1; j++) { 134 for (i = xs; i < xe - 1; i++) { 135 cell[0] = (i - Xs) + (j - Ys) * (Xe - Xs) + (k - Zs) * (Xe - Xs) * (Ye - Ys); 136 cell[1] = (i + 1 - Xs) + (j - Ys) * (Xe - Xs) + (k - Zs) * (Xe - Xs) * (Ye - Ys); 137 cell[2] = (i + 1 - Xs) + (j + 1 - Ys) * (Xe - Xs) + (k - Zs) * (Xe - Xs) * (Ye - Ys); 138 cell[3] = (i - Xs) + (j + 1 - Ys) * (Xe - Xs) + (k - Zs) * (Xe - Xs) * (Ye - Ys); 139 cell[4] = (i - Xs) + (j - Ys) * (Xe - Xs) + (k + 1 - Zs) * (Xe - Xs) * (Ye - Ys); 140 cell[5] = (i + 1 - Xs) + (j - Ys) * (Xe - Xs) + (k + 1 - Zs) * (Xe - Xs) * (Ye - Ys); 141 cell[6] = (i + 1 - Xs) + (j + 1 - Ys) * (Xe - Xs) + (k + 1 - Zs) * (Xe - Xs) * (Ye - Ys); 142 cell[7] = (i - Xs) + (j + 1 - Ys) * (Xe - Xs) + (k + 1 - Zs) * (Xe - Xs) * (Ye - Ys); 143 if (da->elementtype == DMDA_ELEMENT_P1) { 144 for (c = 0; c < ns * nn; c++) da->e[cnt++] = cell[split[c]]; 145 } 146 if (da->elementtype == DMDA_ELEMENT_Q1) { 147 for (c = 0; c < ns * nn; c++) da->e[cnt++] = cell[c]; 148 } 149 } 150 } 151 } 152 153 corners[0] = (xs - Xs) + (ys - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys); 154 corners[1] = (xe - 1 - Xs) + (ys - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys); 155 corners[2] = (xs - Xs) + (ye - 1 - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys); 156 corners[3] = (xe - 1 - Xs) + (ye - 1 - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys); 157 corners[4] = (xs - Xs) + (ys - Ys) * (Xe - Xs) + (ze - 1 - Zs) * (Xe - Xs) * (Ye - Ys); 158 corners[5] = (xe - 1 - Xs) + (ys - Ys) * (Xe - Xs) + (ze - 1 - Zs) * (Xe - Xs) * (Ye - Ys); 159 corners[6] = (xs - Xs) + (ye - 1 - Ys) * (Xe - Xs) + (ze - 1 - Zs) * (Xe - Xs) * (Ye - Ys); 160 corners[7] = (xe - 1 - Xs) + (ye - 1 - Ys) * (Xe - Xs) + (ze - 1 - Zs) * (Xe - Xs) * (Ye - Ys); 161 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 8, corners, PETSC_COPY_VALUES, &da->ecorners)); 162 } 163 *nel = da->ne; 164 *nen = da->nen; 165 *e = da->e; 166 PetscFunctionReturn(0); 167 } 168 169 /*@ 170 DMDAGetElementsCorners - Returns the global (x,y,z) indices of the lower left 171 corner of the non-overlapping decomposition identified by DMDAGetElements() 172 173 Not Collective 174 175 Input Parameter: 176 . da - the DM object 177 178 Output Parameters: 179 + gx - the x index 180 . gy - the y index 181 - gz - the z index 182 183 Level: intermediate 184 185 Notes: 186 187 .seealso: `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()` 188 @*/ 189 PetscErrorCode DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz) { 190 PetscInt xs, Xs; 191 PetscInt ys, Ys; 192 PetscInt zs, Zs; 193 PetscBool isda; 194 195 PetscFunctionBegin; 196 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 197 if (gx) PetscValidIntPointer(gx, 2); 198 if (gy) PetscValidIntPointer(gy, 3); 199 if (gz) PetscValidIntPointer(gz, 4); 200 PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda)); 201 PetscCheck(isda, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)da)->type_name); 202 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, NULL, NULL, NULL)); 203 PetscCall(DMDAGetGhostCorners(da, &Xs, &Ys, &Zs, NULL, NULL, NULL)); 204 if (xs != Xs) xs -= 1; 205 if (ys != Ys) ys -= 1; 206 if (zs != Zs) zs -= 1; 207 if (gx) *gx = xs; 208 if (gy) *gy = ys; 209 if (gz) *gz = zs; 210 PetscFunctionReturn(0); 211 } 212 213 /*@ 214 DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by DMDAGetElements() 215 216 Not Collective 217 218 Input Parameter: 219 . da - the DM object 220 221 Output Parameters: 222 + mx - number of local elements in x-direction 223 . my - number of local elements in y-direction 224 - mz - number of local elements in z-direction 225 226 Level: intermediate 227 228 Notes: 229 It returns the same number of elements, irrespective of the DMDAElementType 230 231 .seealso: `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements` 232 @*/ 233 PetscErrorCode DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz) { 234 PetscInt xs, xe, Xs; 235 PetscInt ys, ye, Ys; 236 PetscInt zs, ze, Zs; 237 PetscInt dim; 238 PetscBool isda; 239 240 PetscFunctionBegin; 241 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 242 if (mx) PetscValidIntPointer(mx, 2); 243 if (my) PetscValidIntPointer(my, 3); 244 if (mz) PetscValidIntPointer(mz, 4); 245 PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda)); 246 PetscCheck(isda, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)da)->type_name); 247 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xe, &ye, &ze)); 248 PetscCall(DMDAGetGhostCorners(da, &Xs, &Ys, &Zs, NULL, NULL, NULL)); 249 xe += xs; 250 if (xs != Xs) xs -= 1; 251 ye += ys; 252 if (ys != Ys) ys -= 1; 253 ze += zs; 254 if (zs != Zs) zs -= 1; 255 if (mx) *mx = 0; 256 if (my) *my = 0; 257 if (mz) *mz = 0; 258 PetscCall(DMGetDimension(da, &dim)); 259 switch (dim) { 260 case 3: 261 if (mz) *mz = ze - zs - 1; 262 case 2: 263 if (my) *my = ye - ys - 1; 264 case 1: 265 if (mx) *mx = xe - xs - 1; 266 break; 267 } 268 PetscFunctionReturn(0); 269 } 270 271 /*@ 272 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 273 274 Not Collective 275 276 Input Parameter: 277 . da - the DMDA object 278 279 Output Parameters: 280 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 281 282 Level: intermediate 283 284 .seealso: `DMDAElementType`, `DMDAGetElementType()`, `DMDAGetElements()`, `DMDARestoreElements()` 285 @*/ 286 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) { 287 DM_DA *dd = (DM_DA *)da->data; 288 PetscBool isda; 289 290 PetscFunctionBegin; 291 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 292 PetscValidLogicalCollectiveEnum(da, etype, 2); 293 PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda)); 294 if (!isda) PetscFunctionReturn(0); 295 if (dd->elementtype != etype) { 296 PetscCall(PetscFree(dd->e)); 297 PetscCall(ISDestroy(&dd->ecorners)); 298 299 dd->elementtype = etype; 300 dd->ne = 0; 301 dd->nen = 0; 302 dd->e = NULL; 303 } 304 PetscFunctionReturn(0); 305 } 306 307 /*@ 308 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 309 310 Not Collective 311 312 Input Parameter: 313 . da - the DMDA object 314 315 Output Parameters: 316 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 317 318 Level: intermediate 319 320 .seealso: `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`, `DMDARestoreElements()` 321 @*/ 322 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) { 323 DM_DA *dd = (DM_DA *)da->data; 324 PetscBool isda; 325 326 PetscFunctionBegin; 327 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 328 PetscValidPointer(etype, 2); 329 PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda)); 330 PetscCheck(isda, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)da)->type_name); 331 *etype = dd->elementtype; 332 PetscFunctionReturn(0); 333 } 334 335 /*@C 336 DMDAGetElements - Gets an array containing the indices (in local coordinates) 337 of all the local elements 338 339 Not Collective 340 341 Input Parameter: 342 . dm - the DM object 343 344 Output Parameters: 345 + nel - number of local elements 346 . nen - number of element nodes 347 - e - the local indices of the elements' vertices 348 349 Level: intermediate 350 351 Notes: 352 Call DMDARestoreElements() once you have finished accessing the elements. 353 354 Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes. 355 356 If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result. 357 358 Not supported in Fortran 359 360 .seealso: `DMDAElementType`, `DMDASetElementType()`, `VecSetValuesLocal()`, `MatSetValuesLocal()`, `DMGlobalToLocalBegin()`, `DMLocalToGlobalBegin()` 361 @*/ 362 PetscErrorCode DMDAGetElements(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[]) { 363 PetscInt dim; 364 DM_DA *dd = (DM_DA *)dm->data; 365 PetscBool isda; 366 367 PetscFunctionBegin; 368 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMDA); 369 PetscValidIntPointer(nel, 2); 370 PetscValidIntPointer(nen, 3); 371 PetscValidPointer(e, 4); 372 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda)); 373 PetscCheck(isda, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)dm)->type_name); 374 PetscCheck(dd->stencil_type != DMDA_STENCIL_STAR, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DMDAGetElements() requires you use a stencil type of DMDA_STENCIL_BOX"); 375 PetscCall(DMGetDimension(dm, &dim)); 376 if (dd->e) { 377 *nel = dd->ne; 378 *nen = dd->nen; 379 *e = dd->e; 380 PetscFunctionReturn(0); 381 } 382 if (dim == -1) { 383 *nel = 0; 384 *nen = 0; 385 *e = NULL; 386 } else if (dim == 1) { 387 PetscCall(DMDAGetElements_1D(dm, nel, nen, e)); 388 } else if (dim == 2) { 389 PetscCall(DMDAGetElements_2D(dm, nel, nen, e)); 390 } else if (dim == 3) { 391 PetscCall(DMDAGetElements_3D(dm, nel, nen, e)); 392 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "DMDA dimension not 1, 2, or 3, it is %" PetscInt_FMT, dim); 393 PetscFunctionReturn(0); 394 } 395 396 /*@ 397 DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates) 398 of the non-overlapping decomposition identified by DMDAGetElements 399 400 Not Collective 401 402 Input Parameter: 403 . dm - the DM object 404 405 Output Parameters: 406 . is - the index set 407 408 Level: intermediate 409 410 Notes: 411 Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set. 412 413 .seealso: `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`, `DMDARestoreElementsCornersIS()` 414 @*/ 415 PetscErrorCode DMDAGetSubdomainCornersIS(DM dm, IS *is) { 416 DM_DA *dd = (DM_DA *)dm->data; 417 PetscBool isda; 418 419 PetscFunctionBegin; 420 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMDA); 421 PetscValidPointer(is, 2); 422 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda)); 423 PetscCheck(isda, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)dm)->type_name); 424 PetscCheck(dd->stencil_type != DMDA_STENCIL_STAR, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DMDAGetElement() requires you use a stencil type of DMDA_STENCIL_BOX"); 425 if (!dd->ecorners) { /* compute elements if not yet done */ 426 const PetscInt *e; 427 PetscInt nel, nen; 428 429 PetscCall(DMDAGetElements(dm, &nel, &nen, &e)); 430 PetscCall(DMDARestoreElements(dm, &nel, &nen, &e)); 431 } 432 *is = dd->ecorners; 433 PetscFunctionReturn(0); 434 } 435 436 /*@C 437 DMDARestoreElements - Restores the array obtained with DMDAGetElements() 438 439 Not Collective 440 441 Input Parameters: 442 + dm - the DM object 443 . nel - number of local elements 444 . nen - number of element nodes 445 - e - the local indices of the elements' vertices 446 447 Level: intermediate 448 449 Note: You should not access these values after you have called this routine. 450 451 This restore signals the DMDA object that you no longer need access to the array information. 452 453 Not supported in Fortran 454 455 .seealso: `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()` 456 @*/ 457 PetscErrorCode DMDARestoreElements(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[]) { 458 PetscFunctionBegin; 459 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMDA); 460 PetscValidIntPointer(nel, 2); 461 PetscValidIntPointer(nen, 3); 462 PetscValidPointer(e, 4); 463 *nel = 0; 464 *nen = -1; 465 *e = NULL; 466 PetscFunctionReturn(0); 467 } 468 469 /*@ 470 DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS() 471 472 Not Collective 473 474 Input Parameters: 475 + dm - the DM object 476 - is - the index set 477 478 Level: intermediate 479 480 Note: 481 482 .seealso: `DMDAElementType`, `DMDASetElementType()`, `DMDAGetSubdomainCornersIS()` 483 @*/ 484 PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm, IS *is) { 485 PetscFunctionBegin; 486 PetscValidHeaderSpecificType(dm, DM_CLASSID, 1, DMDA); 487 PetscValidHeaderSpecific(*is, IS_CLASSID, 2); 488 *is = NULL; 489 PetscFunctionReturn(0); 490 } 491