1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 8 /* 9 Gets the natural number for each global number on the process. 10 11 Used by DMDAGetAO() and DMDAGlobalToNatural_Create() 12 */ 13 PetscErrorCode DMDAGetNatural_Private(DM da, PetscInt *outNlocal, IS *isnatural) { 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(0); 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 Notes: 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: `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 DM_DA *dd; 70 PetscBool isdmda; 71 72 PetscFunctionBegin; 73 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 74 PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isdmda)); 75 PetscCheck(isdmda, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Requires a DMDA as input"); 76 /* now we can safely dereference */ 77 dd = (DM_DA *)da->data; 78 if (dd->ao) { /* check if the already computed AO has the same type as requested */ 79 PetscBool match; 80 PetscCall(PetscObjectTypeCompare((PetscObject)dd->ao, aotype, &match)); 81 PetscCheck(match, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Cannot change AO type"); 82 PetscFunctionReturn(0); 83 } 84 PetscCall(PetscFree(dd->aotype)); 85 PetscCall(PetscStrallocpy(aotype, (char **)&dd->aotype)); 86 PetscFunctionReturn(0); 87 } 88 89 /*@ 90 DMDAGetAO - Gets the application ordering context for a distributed array. 91 92 Collective on da 93 94 Input Parameter: 95 . da - the distributed array 96 97 Output Parameters: 98 . ao - the application ordering context for DMDAs 99 100 Level: intermediate 101 102 Notes: 103 In this case, the AO maps to the natural grid ordering that would be used 104 for the DMDA if only 1 processor were employed (ordering most rapidly in the 105 x-direction, then y, then z). Multiple degrees of freedom are numbered 106 for each node (rather than 1 component for the whole grid, then the next 107 component, etc.) 108 109 Do NOT call AODestroy() on the ao returned by this function. 110 111 .seealso: `DMDACreate2d()`, `DMDASetAOType()`, `DMDAGetGhostCorners()`, `DMDAGetCorners()`, `DMLocalToGlobal()` 112 `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDAGetOwnershipRanges()`, 113 `AO`, `AOPetscToApplication()`, `AOApplicationToPetsc()` 114 @*/ 115 PetscErrorCode DMDAGetAO(DM da, AO *ao) { 116 DM_DA *dd; 117 PetscBool isdmda; 118 119 PetscFunctionBegin; 120 PetscValidHeaderSpecificType(da, DM_CLASSID, 1, DMDA); 121 PetscValidPointer(ao, 2); 122 PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isdmda)); 123 PetscCheck(isdmda, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Requires a DMDA as input"); 124 /* now we can safely dereference */ 125 dd = (DM_DA *)da->data; 126 127 /* 128 Build the natural ordering to PETSc ordering mappings. 129 */ 130 if (!dd->ao) { 131 IS ispetsc, isnatural; 132 PetscInt Nlocal; 133 134 PetscCall(DMDAGetNatural_Private(da, &Nlocal, &isnatural)); 135 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)da), Nlocal, dd->base, 1, &ispetsc)); 136 PetscCall(AOCreate(PetscObjectComm((PetscObject)da), &dd->ao)); 137 PetscCall(AOSetIS(dd->ao, isnatural, ispetsc)); 138 PetscCall(AOSetType(dd->ao, dd->aotype)); 139 PetscCall(ISDestroy(&ispetsc)); 140 PetscCall(ISDestroy(&isnatural)); 141 } 142 *ao = dd->ao; 143 PetscFunctionReturn(0); 144 } 145