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