1: /*
  2:   Code for manipulating distributed regular arrays in parallel.
  3: */

  5: #include <petsc/private/dmdaimpl.h>

  7: /*
  8:    Gets the natural number for each global number on the process.

 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;

 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);

 22:   PetscCall(PetscMalloc1(Nlocal, &lidx));

 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: }

 48: /*@
 49:   DMDASetAOType - Sets the type of application ordering to create with `DMDAGetAO()`, for a distributed array.

 51:   Collective

 53:   Input Parameters:
 54: + da     - the `DMDA`
 55: - aotype - type of `AO`. `AOType` which can be `AOBASIC`, `AOADVANCED`, `AOMAPPING`, or `AOMEMORYSCALABLE`

 57:   Level: intermediate

 59:   Note:
 60:   It will generate an error if an `AO` has already been obtained with a call to `DMDAGetAO()` and the user sets a different `AOType`

 62: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate2d()`, `DMDAGetAO()`, `DMDAGetGhostCorners()`, `DMDAGetCorners()`, `DMLocalToGlobal()`
 63:           `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDAGetGlobalIndices()`, `DMDAGetOwnershipRanges()`,
 64:           `AO`, `AOPetscToApplication()`, `AOApplicationToPetsc()`, `AOType`, `AOBASIC`, `AOADVANCED`, `AOMAPPING`, `AOMEMORYSCALABLE`
 65: @*/
 66: PetscErrorCode DMDASetAOType(DM da, AOType aotype)
 67: {
 68:   DM_DA    *dd;
 69:   PetscBool isdmda;

 71:   PetscFunctionBegin;
 73:   PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isdmda));
 74:   PetscCheck(isdmda, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Requires a DMDA as input");
 75:   /* now we can safely dereference */
 76:   dd = (DM_DA *)da->data;
 77:   if (dd->ao) { /* check if the already computed AO has the same type as requested */
 78:     PetscBool match;
 79:     PetscCall(PetscObjectTypeCompare((PetscObject)dd->ao, aotype, &match));
 80:     PetscCheck(match, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Cannot change AO type");
 81:     PetscFunctionReturn(PETSC_SUCCESS);
 82:   }
 83:   PetscCall(PetscFree(dd->aotype));
 84:   PetscCall(PetscStrallocpy(aotype, (char **)&dd->aotype));
 85:   PetscFunctionReturn(PETSC_SUCCESS);
 86: }

 88: /*@
 89:   DMDAGetAO - Gets the application ordering context for a distributed array.

 91:   Collective

 93:   Input Parameter:
 94: . da - the `DMDA`

 96:   Output Parameter:
 97: . ao - the application ordering context for `DMDA`

 99:   Level: intermediate

101:   Notes:
102:   In this case, the `AO` maps to the natural grid ordering that would be used
103:   for the `DMDA` if only 1 processor were employed (ordering most rapidly in the
104:   x-direction, then y, then z).  Multiple degrees of freedom are numbered
105:   for each node (rather than 1 component for the whole grid, then the next
106:   component, etc.)

108:   Do NOT call `AODestroy()` on the `ao` returned by this function.

110: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate2d()`, `DMDASetAOType()`, `DMDAGetGhostCorners()`, `DMDAGetCorners()`, `DMLocalToGlobal()`
111:           `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDAGetOwnershipRanges()`,
112:           `AO`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
113: @*/
114: PetscErrorCode DMDAGetAO(DM da, AO *ao)
115: {
116:   DM_DA    *dd;
117:   PetscBool isdmda;

119:   PetscFunctionBegin;
121:   PetscAssertPointer(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;

127:   /*
128:      Build the natural ordering to PETSc ordering mappings.
129:   */
130:   if (!dd->ao) {
131:     IS       ispetsc, isnatural;
132:     PetscInt Nlocal;

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(PETSC_SUCCESS);
144: }