Actual source code: daindex.c

  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:   Nlocal = (dd->xe - dd->xs);
 18:   if (dim > 1) Nlocal *= (dd->ye - dd->ys);
 19:   if (dim > 2) Nlocal *= (dd->ze - dd->zs);

 21:   PetscMalloc1(Nlocal, &lidx);

 23:   if (dim == 1) {
 24:     for (i = dd->xs; i < dd->xe; i++) {
 25:       /*  global number in natural ordering */
 26:       lidx[lict++] = i;
 27:     }
 28:   } else if (dim == 2) {
 29:     for (j = dd->ys; j < dd->ye; j++) {
 30:       for (i = dd->xs; i < dd->xe; i++) {
 31:         /*  global number in natural ordering */
 32:         lidx[lict++] = i + j * dd->M * dd->w;
 33:       }
 34:     }
 35:   } else if (dim == 3) {
 36:     for (k = dd->zs; k < dd->ze; k++) {
 37:       for (j = dd->ys; j < dd->ye; j++) {
 38:         for (i = dd->xs; i < dd->xe; i++) lidx[lict++] = i + j * dd->M * dd->w + k * dd->M * dd->N * dd->w;
 39:       }
 40:     }
 41:   }
 42:   *outNlocal = Nlocal;
 43:   ISCreateGeneral(PetscObjectComm((PetscObject)da), Nlocal, lidx, PETSC_OWN_POINTER, isnatural);
 44:   return 0;
 45: }

 47: /*@C
 48:    DMDASetAOType - Sets the type of application ordering for a distributed array.

 50:    Collective on da

 52:    Input Parameters:
 53: +  da - the distributed array
 54: -  aotype - type of `AO`

 56:    Output Parameters:

 58:    Level: intermediate

 60:    Note:
 61:    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: .seealso: `DM`, `DMDA`, `DMDACreate2d()`, `DMDAGetAO()`, `DMDAGetGhostCorners()`, `DMDAGetCorners()`, `DMLocalToGlobal()`
 64:           `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDAGetGlobalIndices()`, `DMDAGetOwnershipRanges()`,
 65:           `AO`, `AOPetscToApplication()`, `AOApplicationToPetsc()`
 66: @*/
 67: PetscErrorCode DMDASetAOType(DM da, AOType aotype)
 68: {
 69:   DM_DA    *dd;
 70:   PetscBool isdmda;

 73:   PetscObjectTypeCompare((PetscObject)da, DMDA, &isdmda);
 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:     PetscObjectTypeCompare((PetscObject)dd->ao, aotype, &match);
 81:     return 0;
 82:   }
 83:   PetscFree(dd->aotype);
 84:   PetscStrallocpy(aotype, (char **)&dd->aotype);
 85:   return 0;
 86: }

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

 91:    Collective on da

 93:    Input Parameter:
 94: .  da - the distributed array

 96:    Output Parameters:
 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: `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;

121:   PetscObjectTypeCompare((PetscObject)da, DMDA, &isdmda);
123:   /* now we can safely dereference */
124:   dd = (DM_DA *)da->data;

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

133:     DMDAGetNatural_Private(da, &Nlocal, &isnatural);
134:     ISCreateStride(PetscObjectComm((PetscObject)da), Nlocal, dd->base, 1, &ispetsc);
135:     AOCreate(PetscObjectComm((PetscObject)da), &dd->ao);
136:     AOSetIS(dd->ao, isnatural, ispetsc);
137:     AOSetType(dd->ao, dd->aotype);
138:     ISDestroy(&ispetsc);
139:     ISDestroy(&isnatural);
140:   }
141:   *ao = dd->ao;
142:   return 0;
143: }