1 /* 2 Code for manipulating distributed regular arrays in parallel. 3 */ 4 5 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 6 7 /* 8 Gets the natural number for each global number on the process. 9 10 Used by DMDAGetAO() and DMDAGlobalToNatural_Create() 11 */ 12 PetscErrorCode DMDAGetNatural_Private(DM da, PetscInt *outNlocal, IS *isnatural) 13 { 14 PetscInt Nlocal, i, j, k, *lidx, lict = 0, dim = da->dim; 15 DM_DA *dd = (DM_DA *)da->data; 16 17 PetscFunctionBegin; 18 Nlocal = (dd->xe - dd->xs); 19 if (dim > 1) Nlocal *= (dd->ye - dd->ys); 20 if (dim > 2) Nlocal *= (dd->ze - dd->zs); 21 22 PetscCall(PetscMalloc1(Nlocal, &lidx)); 23 24 if (dim == 1) { 25 for (i = dd->xs; i < dd->xe; i++) { 26 /* global number in natural ordering */ 27 lidx[lict++] = i; 28 } 29 } else if (dim == 2) { 30 for (j = dd->ys; j < dd->ye; j++) { 31 for (i = dd->xs; i < dd->xe; i++) { 32 /* global number in natural ordering */ 33 lidx[lict++] = i + j * dd->M * dd->w; 34 } 35 } 36 } else if (dim == 3) { 37 for (k = dd->zs; k < dd->ze; k++) { 38 for (j = dd->ys; j < dd->ye; j++) { 39 for (i = dd->xs; i < dd->xe; i++) lidx[lict++] = i + j * dd->M * dd->w + k * dd->M * dd->N * dd->w; 40 } 41 } 42 } 43 *outNlocal = Nlocal; 44 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)da), Nlocal, lidx, PETSC_OWN_POINTER, isnatural)); 45 PetscFunctionReturn(PETSC_SUCCESS); 46 } 47 48 /*@C 49 DMDASetAOType - Sets the type of application ordering for a distributed array. 50 51 Collective on da 52 53 Input Parameters: 54 + da - the distributed array 55 - aotype - type of `AO` 56 57 Output Parameters: 58 59 Level: intermediate 60 61 Note: 62 It will generate and error if an `AO` has already been obtained with a call to `DMDAGetAO()` and the user sets a different `AOType` 63 64 .seealso: `DM`, `DMDA`, `DMDACreate2d()`, `DMDAGetAO()`, `DMDAGetGhostCorners()`, `DMDAGetCorners()`, `DMLocalToGlobal()` 65 `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDAGetGlobalIndices()`, `DMDAGetOwnershipRanges()`, 66 `AO`, `AOPetscToApplication()`, `AOApplicationToPetsc()` 67 @*/ 68 PetscErrorCode DMDASetAOType(DM da, AOType aotype) 69 { 70 DM_DA *dd; 71 PetscBool isdmda; 72 73 PetscFunctionBegin; 74 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 75 PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isdmda)); 76 PetscCheck(isdmda, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Requires a DMDA as input"); 77 /* now we can safely dereference */ 78 dd = (DM_DA *)da->data; 79 if (dd->ao) { /* check if the already computed AO has the same type as requested */ 80 PetscBool match; 81 PetscCall(PetscObjectTypeCompare((PetscObject)dd->ao, aotype, &match)); 82 PetscCheck(match, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Cannot change AO type"); 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 PetscCall(PetscFree(dd->aotype)); 86 PetscCall(PetscStrallocpy(aotype, (char **)&dd->aotype)); 87 PetscFunctionReturn(PETSC_SUCCESS); 88 } 89 90 /*@ 91 DMDAGetAO - Gets the application ordering context for a distributed array. 92 93 Collective on da 94 95 Input Parameter: 96 . da - the distributed array 97 98 Output Parameters: 99 . ao - the application ordering context for `DMDA` 100 101 Level: intermediate 102 103 Notes: 104 In this case, the `AO` maps to the natural grid ordering that would be used 105 for the `DMDA` if only 1 processor were employed (ordering most rapidly in the 106 x-direction, then y, then z). Multiple degrees of freedom are numbered 107 for each node (rather than 1 component for the whole grid, then the next 108 component, etc.) 109 110 Do NOT call `AODestroy()` on the ao returned by this function. 111 112 .seealso: `DM`, `DMDA`, `DMDACreate2d()`, `DMDASetAOType()`, `DMDAGetGhostCorners()`, `DMDAGetCorners()`, `DMLocalToGlobal()` 113 `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDAGetOwnershipRanges()`, 114 `AO`, `AOPetscToApplication()`, `AOApplicationToPetsc()` 115 @*/ 116 PetscErrorCode DMDAGetAO(DM da, AO *ao) 117 { 118 DM_DA *dd; 119 PetscBool isdmda; 120 121 PetscFunctionBegin; 122 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 123 PetscValidPointer(ao, 2); 124 PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isdmda)); 125 PetscCheck(isdmda, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Requires a DMDA as input"); 126 /* now we can safely dereference */ 127 dd = (DM_DA *)da->data; 128 129 /* 130 Build the natural ordering to PETSc ordering mappings. 131 */ 132 if (!dd->ao) { 133 IS ispetsc, isnatural; 134 PetscInt Nlocal; 135 136 PetscCall(DMDAGetNatural_Private(da, &Nlocal, &isnatural)); 137 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)da), Nlocal, dd->base, 1, &ispetsc)); 138 PetscCall(AOCreate(PetscObjectComm((PetscObject)da), &dd->ao)); 139 PetscCall(AOSetIS(dd->ao, isnatural, ispetsc)); 140 PetscCall(AOSetType(dd->ao, dd->aotype)); 141 PetscCall(ISDestroy(&ispetsc)); 142 PetscCall(ISDestroy(&isnatural)); 143 } 144 *ao = dd->ao; 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147