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