Actual source code: daindex.c


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

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

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

 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;

 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:   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:   ISCreateGeneral(PetscObjectComm((PetscObject)da), Nlocal, lidx, PETSC_OWN_POINTER, isnatural);
 45:   return 0;
 46: }

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

 51:    Collective on da

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

 57:    Output Parameters:

 59:    Level: intermediate

 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

 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: {
 70:   DM_DA    *dd;
 71:   PetscBool isdmda;

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

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

 92:    Collective on da

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

 97:    Output Parameters:
 98: .  ao - the application ordering context for DMDAs

100:    Level: intermediate

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

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

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: {
117:   DM_DA    *dd;
118:   PetscBool isdmda;

122:   PetscObjectTypeCompare((PetscObject)da, DMDA, &isdmda);
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:     DMDAGetNatural_Private(da, &Nlocal, &isnatural);
135:     ISCreateStride(PetscObjectComm((PetscObject)da), Nlocal, dd->base, 1, &ispetsc);
136:     AOCreate(PetscObjectComm((PetscObject)da), &dd->ao);
137:     AOSetIS(dd->ao, isnatural, ispetsc);
138:     AOSetType(dd->ao, dd->aotype);
139:     ISDestroy(&ispetsc);
140:     ISDestroy(&isnatural);
141:   }
142:   *ao = dd->ao;
143:   return 0;
144: }