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