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: {
 16:   PetscInt       Nlocal,i,j,k,*lidx,lict = 0,dim = da->dim;
 17:   DM_DA          *dd = (DM_DA*)da->data;

 20:   Nlocal = (dd->xe-dd->xs);
 21:   if (dim > 1) Nlocal *= (dd->ye-dd->ys);
 22:   if (dim > 2) Nlocal *= (dd->ze-dd->zs);

 24:   PetscMalloc1(Nlocal,&lidx);

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

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

 55:    Collective on da

 57:    Input Parameters:
 58: +  da - the distributed array
 59: -  aotype - type of AO

 61:    Output Parameters:

 63:    Level: intermediate

 65:    Notes:
 66:    It will generate and error if an AO has already been obtained with a call to DMDAGetAO and the user sets a different AOType

 68: .seealso: DMDACreate2d(), DMDAGetAO(), DMDAGetGhostCorners(), DMDAGetCorners(), DMLocalToGlobal()
 69:           DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDAGetGlobalIndices(), DMDAGetOwnershipRanges(),
 70:           AO, AOPetscToApplication(), AOApplicationToPetsc()
 71: @*/
 72: PetscErrorCode  DMDASetAOType(DM da,AOType aotype)
 73: {
 74:   DM_DA          *dd;
 75:   PetscBool      isdmda;

 80:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isdmda);
 81:   if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Requires a DMDA as input");
 82:   /* now we can safely dereference */
 83:   dd = (DM_DA*)da->data;
 84:   if (dd->ao) { /* check if the already computed AO has the same type as requested */
 85:     PetscBool match;
 86:     PetscObjectTypeCompare((PetscObject)dd->ao,aotype,&match);
 87:     if (!match) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot change AO type");
 88:     return(0);
 89:   }
 90:   PetscFree(dd->aotype);
 91:   PetscStrallocpy(aotype,(char**)&dd->aotype);
 92:   return(0);
 93: }

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

 98:    Collective on da

100:    Input Parameter:
101: .  da - the distributed array

103:    Output Parameters:
104: .  ao - the application ordering context for DMDAs

106:    Level: intermediate

108:    Notes:
109:    In this case, the AO maps to the natural grid ordering that would be used
110:    for the DMDA if only 1 processor were employed (ordering most rapidly in the
111:    x-direction, then y, then z).  Multiple degrees of freedom are numbered
112:    for each node (rather than 1 component for the whole grid, then the next
113:    component, etc.)

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

117: .seealso: DMDACreate2d(), DMDASetAOType(), DMDAGetGhostCorners(), DMDAGetCorners(), DMLocalToGlobal()
118:           DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToLocalBegin(), DMLocalToLocalEnd(),  DMDAGetOwnershipRanges(),
119:           AO, AOPetscToApplication(), AOApplicationToPetsc()
120: @*/
121: PetscErrorCode  DMDAGetAO(DM da,AO *ao)
122: {
123:   DM_DA          *dd;
124:   PetscBool      isdmda;

130:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isdmda);
131:   if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Requires a DMDA as input");
132:   /* now we can safely dereference */
133:   dd = (DM_DA*)da->data;

135:   /*
136:      Build the natural ordering to PETSc ordering mappings.
137:   */
138:   if (!dd->ao) {
139:     IS             ispetsc,isnatural;
141:     PetscInt       Nlocal;

143:     DMDAGetNatural_Private(da,&Nlocal,&isnatural);
144:     ISCreateStride(PetscObjectComm((PetscObject)da),Nlocal,dd->base,1,&ispetsc);
145:     AOCreate(PetscObjectComm((PetscObject)da),&dd->ao);
146:     AOSetIS(dd->ao,isnatural,ispetsc);
147:     AOSetType(dd->ao,dd->aotype);
148:     PetscLogObjectParent((PetscObject)da,(PetscObject)dd->ao);
149:     ISDestroy(&ispetsc);
150:     ISDestroy(&isnatural);
151:   }
152:   *ao = dd->ao;
153:   return(0);
154: }