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

  6: #include <petsc/private/dmdaimpl.h>    /*I   "petscdmda.h"   I*/

 10: /*
 11:    Gets the natural number for each global number on the process.

 13:    Used by DMDAGetAO() and DMDAGlobalToNatural_Create()
 14: */
 15: PetscErrorCode DMDAGetNatural_Private(DM da,PetscInt *outNlocal,IS *isnatural)
 16: {
 18:   PetscInt       Nlocal,i,j,k,*lidx,lict = 0,dim = da->dim;
 19:   DM_DA          *dd = (DM_DA*)da->data;

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

 26:   PetscMalloc1(Nlocal,&lidx);

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

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

 59:    Collective on DMDA

 61:    Input Parameter:
 62: .  da - the distributed array
 63: .  aotype - type of AO

 65:    Output Parameters:

 67:    Level: intermediate

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

 72: .keywords: distributed array, get, global, indices, local-to-global

 74: .seealso: DMDACreate2d(), DMDAGetAO(), DMDAGetGhostCorners(), DMDAGetCorners(), DMDALocalToGlocal()
 75:           DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDAGetGlobalIndices(), DMDAGetOwnershipRanges(),
 76:           AO, AOPetscToApplication(), AOApplicationToPetsc()
 77: @*/
 78: PetscErrorCode  DMDASetAOType(DM da,AOType aotype)
 79: {
 80:   DM_DA          *dd;
 81:   PetscBool      isdmda;

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

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

106:    Collective on DMDA

108:    Input Parameter:
109: .  da - the distributed array

111:    Output Parameters:
112: .  ao - the application ordering context for DMDAs

114:    Level: intermediate

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

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

125: .keywords: distributed array, get, global, indices, local-to-global

127: .seealso: DMDACreate2d(), DMDASetAOType(), DMDAGetGhostCorners(), DMDAGetCorners(), DMDALocalToGlocal()
128:           DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToLocalBegin(), DMLocalToLocalEnd(),  DMDAGetOwnershipRanges(),
129:           AO, AOPetscToApplication(), AOApplicationToPetsc()
130: @*/
131: PetscErrorCode  DMDAGetAO(DM da,AO *ao)
132: {
133:   DM_DA          *dd;
134:   PetscBool      isdmda;

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

145:   /*
146:      Build the natural ordering to PETSc ordering mappings.
147:   */
148:   if (!dd->ao) {
149:     IS             ispetsc,isnatural;
151:     PetscInt       Nlocal;

153:     DMDAGetNatural_Private(da,&Nlocal,&isnatural);
154:     ISCreateStride(PetscObjectComm((PetscObject)da),Nlocal,dd->base,1,&ispetsc);
155:     AOCreate(PetscObjectComm((PetscObject)da),&dd->ao);
156:     AOSetIS(dd->ao,isnatural,ispetsc);
157:     AOSetType(dd->ao,dd->aotype);
158:     PetscLogObjectParent((PetscObject)da,(PetscObject)dd->ao);
159:     ISDestroy(&ispetsc);
160:     ISDestroy(&isnatural);
161:   }
162:   *ao = dd->ao;
163:   return(0);
164: }