Actual source code: da.c

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

      3: /*@
      4:   DMDASetSizes - Sets the number of grid points in the three dimensional directions

      6:   Logically Collective

      8:   Input Parameters:
      9: + da - the `DMDA`
     10: . M  - the global X size
     11: . N  - the global Y size
     12: - P  - the global Z size

     14:   Level: intermediate

     16:   Developer Note:
     17:   Since the dimension may not yet have been set the code cannot error check for non-positive Y and Z number of grid points

     19: .seealso: [](sec_struct), `DM`, `DMDA`, `PetscSplitOwnership()`
     20: @*/
     21: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
     22: {
     23:   DM_DA *dd = (DM_DA *)da->data;

     25:   PetscFunctionBegin;
     30:   PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
     31:   PetscCheck(M >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Number of grid points in X direction must be positive");
     32:   PetscCheck(N >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Number of grid points in Y direction must be positive");
     33:   PetscCheck(P >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Number of grid points in Z direction must be positive");

     35:   dd->M = M;
     36:   dd->N = N;
     37:   dd->P = P;
     38:   PetscFunctionReturn(PETSC_SUCCESS);
     39: }

     41: /*@
     42:   DMDASetNumProcs - Sets the number of processes in each dimension

     44:   Logically Collective

     46:   Input Parameters:
     47: + da - the `DMDA`
     48: . m  - the number of X processes (or `PETSC_DECIDE`)
     49: . n  - the number of Y processes (or `PETSC_DECIDE`)
     50: - p  - the number of Z processes (or `PETSC_DECIDE`)

     52:   Level: intermediate

     54: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDASetSizes()`, `PetscSplitOwnership()`
     55: @*/
     56: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
     57: {
     58:   DM_DA *dd = (DM_DA *)da->data;

     60:   PetscFunctionBegin;
     65:   PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
     66:   dd->m = m;
     67:   dd->n = n;
     68:   dd->p = p;
     69:   if (da->dim == 2) {
     70:     PetscMPIInt size;
     71:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
     72:     if ((dd->m > 0) && (dd->n < 0)) {
     73:       dd->n = size / dd->m;
     74:       PetscCheck(dd->n * dd->m == size, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " processes in X direction not divisible into comm size %d", m, size);
     75:     }
     76:     if ((dd->n > 0) && (dd->m < 0)) {
     77:       dd->m = size / dd->n;
     78:       PetscCheck(dd->n * dd->m == size, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " processes in Y direction not divisible into comm size %d", n, size);
     79:     }
     80:   }
     81:   PetscFunctionReturn(PETSC_SUCCESS);
     82: }

     84: /*@
     85:   DMDAGetBoundaryType - Gets the type of ghost nodes on domain boundaries.

     87:   Not Collective

     89:   Input Parameter:
     90: . da - The `DMDA`

     92:   Output Parameters:
     93: + bx - x boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
     94: . by - y boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
     95: - bz - z boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`

     97:   Level: intermediate

     99: .seealso: [](sec_struct), `DMDASetBoundaryType()`, `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMBoundaryType`, `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
    100: @*/
    101: PetscErrorCode DMDAGetBoundaryType(DM da, PeOp DMBoundaryType *bx, PeOp DMBoundaryType *by, PeOp DMBoundaryType *bz)
    102: {
    103:   DM_DA *dd = (DM_DA *)da->data;

    105:   PetscFunctionBegin;
    107:   if (bx) PetscAssertPointer(bx, 2);
    108:   if (by) PetscAssertPointer(by, 3);
    109:   if (bz) PetscAssertPointer(bz, 4);
    110:   if (bx) *bx = dd->bx;
    111:   if (by) *by = dd->by;
    112:   if (bz) *bz = dd->bz;
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    114: }

    116: /*@
    117:   DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries for a `DMDA` object.

    119:   Not Collective

    121:   Input Parameters:
    122: + da - The `DMDA`
    123: . bx - x boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
    124: . by - y boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
    125: - bz - z boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`

    127:   Level: intermediate

    129:   Note:
    130:   The default is `DM_BOUNDARY_NONE`

    132: .seealso: [](sec_struct), `DMDAGetBoundaryType()`, `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMBoundaryType`, `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
    133: @*/
    134: PetscErrorCode DMDASetBoundaryType(DM da, DMBoundaryType bx, DMBoundaryType by, DMBoundaryType bz)
    135: {
    136:   DM_DA *dd = (DM_DA *)da->data;

    138:   PetscFunctionBegin;
    143:   PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    144:   dd->bx = bx;
    145:   dd->by = by;
    146:   dd->bz = bz;
    147:   PetscFunctionReturn(PETSC_SUCCESS);
    148: }

    150: /*@
    151:   DMDASetDof - Sets the number of degrees of freedom per vertex

    153:   Not Collective

    155:   Input Parameters:
    156: + da  - The `DMDA`
    157: - dof - Number of degrees of freedom per vertex

    159:   Level: intermediate

    161: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetDof()`, `DMDACreate()`, `DMDestroy()`
    162: @*/
    163: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
    164: {
    165:   DM_DA *dd = (DM_DA *)da->data;

    167:   PetscFunctionBegin;
    170:   PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    171:   dd->w  = dof;
    172:   da->bs = dof;
    173:   PetscFunctionReturn(PETSC_SUCCESS);
    174: }

    176: /*@
    177:   DMDAGetDof - Gets the number of degrees of freedom per vertex

    179:   Not Collective

    181:   Input Parameter:
    182: . da - The `DMDA`

    184:   Output Parameter:
    185: . dof - Number of degrees of freedom per vertex

    187:   Level: intermediate

    189: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDASetDof()`, `DMDACreate()`, `DMDestroy()`
    190: @*/
    191: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
    192: {
    193:   DM_DA *dd = (DM_DA *)da->data;

    195:   PetscFunctionBegin;
    197:   PetscAssertPointer(dof, 2);
    198:   *dof = dd->w;
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: /*@
    203:   DMDAGetOverlap - Gets the size of the per-processor overlap.

    205:   Not Collective

    207:   Input Parameter:
    208: . da - The `DMDA`

    210:   Output Parameters:
    211: + x - Overlap in the x direction
    212: . y - Overlap in the y direction
    213: - z - Overlap in the z direction

    215:   Level: intermediate

    217: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDASetOverlap()`
    218: @*/
    219: PetscErrorCode DMDAGetOverlap(DM da, PeOp PetscInt *x, PeOp PetscInt *y, PeOp PetscInt *z)
    220: {
    221:   DM_DA *dd = (DM_DA *)da->data;

    223:   PetscFunctionBegin;
    225:   if (x) *x = dd->xol;
    226:   if (y) *y = dd->yol;
    227:   if (z) *z = dd->zol;
    228:   PetscFunctionReturn(PETSC_SUCCESS);
    229: }

    231: /*@
    232:   DMDASetOverlap - Sets the size of the per-processor overlap.

    234:   Not Collective

    236:   Input Parameters:
    237: + da - The `DMDA`
    238: . x  - Overlap in the x direction
    239: . y  - Overlap in the y direction
    240: - z  - Overlap in the z direction

    242:   Level: intermediate

    244: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDAGetOverlap()`
    245: @*/
    246: PetscErrorCode DMDASetOverlap(DM da, PetscInt x, PetscInt y, PetscInt z)
    247: {
    248:   DM_DA *dd = (DM_DA *)da->data;

    250:   PetscFunctionBegin;
    255:   dd->xol = x;
    256:   dd->yol = y;
    257:   dd->zol = z;
    258:   PetscFunctionReturn(PETSC_SUCCESS);
    259: }

    261: /*@
    262:   DMDAGetNumLocalSubDomains - Gets the number of local subdomains that would be created upon decomposition.

    264:   Not Collective

    266:   Input Parameter:
    267: . da - The `DMDA`

    269:   Output Parameter:
    270: . Nsub - Number of local subdomains created upon decomposition

    272:   Level: intermediate

    274: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDASetNumLocalSubDomains()`
    275: @*/
    276: PetscErrorCode DMDAGetNumLocalSubDomains(DM da, PetscInt *Nsub)
    277: {
    278:   DM_DA *dd = (DM_DA *)da->data;

    280:   PetscFunctionBegin;
    282:   if (Nsub) *Nsub = dd->Nsub;
    283:   PetscFunctionReturn(PETSC_SUCCESS);
    284: }

    286: /*@
    287:   DMDASetNumLocalSubDomains - Sets the number of local subdomains to create when decomposing with `DMCreateDomainDecomposition()`

    289:   Not Collective

    291:   Input Parameters:
    292: + da   - The `DMDA`
    293: - Nsub - The number of local subdomains requested

    295:   Level: intermediate

    297: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDAGetNumLocalSubDomains()`
    298: @*/
    299: PetscErrorCode DMDASetNumLocalSubDomains(DM da, PetscInt Nsub)
    300: {
    301:   DM_DA *dd = (DM_DA *)da->data;

    303:   PetscFunctionBegin;
    306:   dd->Nsub = Nsub;
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }

    310: /*@
    311:   DMDASetOffset - Sets the index offset of the `DMDA`.

    313:   Collective

    315:   Input Parameters:
    316: + da - The `DMDA`
    317: . xo - The offset in the x direction
    318: . yo - The offset in the y direction
    319: . zo - The offset in the z direction
    320: . Mo - The problem offset in the x direction
    321: . No - The problem offset in the y direction
    322: - Po - The problem offset in the z direction

    324:   Level: developer

    326:   Note:
    327:   This is used primarily to overlap a computation on a local `DMDA` with that on a global `DMDA` without
    328:   changing boundary conditions or subdomain features that depend upon the global offsets.

    330: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
    331: @*/
    332: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
    333: {
    334:   DM_DA *dd = (DM_DA *)da->data;

    336:   PetscFunctionBegin;
    344:   dd->xo = xo;
    345:   dd->yo = yo;
    346:   dd->zo = zo;
    347:   dd->Mo = Mo;
    348:   dd->No = No;
    349:   dd->Po = Po;

    351:   if (da->coordinates[0].dm) PetscCall(DMDASetOffset(da->coordinates[0].dm, xo, yo, zo, Mo, No, Po));
    352:   PetscFunctionReturn(PETSC_SUCCESS);
    353: }

    355: /*@
    356:   DMDAGetOffset - Gets the index offset of the `DMDA`.

    358:   Not Collective

    360:   Input Parameter:
    361: . da - The `DMDA`

    363:   Output Parameters:
    364: + xo - The offset in the x direction
    365: . yo - The offset in the y direction
    366: . zo - The offset in the z direction
    367: . Mo - The global size in the x direction
    368: . No - The global size in the y direction
    369: - Po - The global size in the z direction

    371:   Level: developer

    373: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDASetOffset()`, `DMDAVecGetArray()`
    374: @*/
    375: PetscErrorCode DMDAGetOffset(DM da, PeOp PetscInt *xo, PeOp PetscInt *yo, PeOp PetscInt *zo, PeOp PetscInt *Mo, PeOp PetscInt *No, PeOp PetscInt *Po)
    376: {
    377:   DM_DA *dd = (DM_DA *)da->data;

    379:   PetscFunctionBegin;
    381:   if (xo) *xo = dd->xo;
    382:   if (yo) *yo = dd->yo;
    383:   if (zo) *zo = dd->zo;
    384:   if (Mo) *Mo = dd->Mo;
    385:   if (No) *No = dd->No;
    386:   if (Po) *Po = dd->Po;
    387:   PetscFunctionReturn(PETSC_SUCCESS);
    388: }

    390: /*@
    391:   DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain `DMDA`.

    393:   Not Collective

    395:   Input Parameter:
    396: . da - The `DMDA`

    398:   Output Parameters:
    399: + xs - The start of the region in x
    400: . ys - The start of the region in y
    401: . zs - The start of the region in z
    402: . xm - The size of the region in x
    403: . ym - The size of the region in y
    404: - zm - The size of the region in z

    406:   Level: intermediate

    408: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
    409: @*/
    410: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PeOp PetscInt *xs, PeOp PetscInt *ys, PeOp PetscInt *zs, PeOp PetscInt *xm, PeOp PetscInt *ym, PeOp PetscInt *zm)
    411: {
    412:   DM_DA *dd = (DM_DA *)da->data;

    414:   PetscFunctionBegin;
    416:   if (xs) *xs = dd->nonxs;
    417:   if (ys) *ys = dd->nonys;
    418:   if (zs) *zs = dd->nonzs;
    419:   if (xm) *xm = dd->nonxm;
    420:   if (ym) *ym = dd->nonym;
    421:   if (zm) *zm = dd->nonzm;
    422:   PetscFunctionReturn(PETSC_SUCCESS);
    423: }

    425: /*@
    426:   DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain `DMDA`.

    428:   Collective

    430:   Input Parameters:
    431: + da - The `DMDA`
    432: . xs - The start of the region in x
    433: . ys - The start of the region in y
    434: . zs - The start of the region in z
    435: . xm - The size of the region in x
    436: . ym - The size of the region in y
    437: - zm - The size of the region in z

    439:   Level: intermediate

    441: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
    442: @*/
    443: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
    444: {
    445:   DM_DA *dd = (DM_DA *)da->data;

    447:   PetscFunctionBegin;
    455:   dd->nonxs = xs;
    456:   dd->nonys = ys;
    457:   dd->nonzs = zs;
    458:   dd->nonxm = xm;
    459:   dd->nonym = ym;
    460:   dd->nonzm = zm;
    461:   PetscFunctionReturn(PETSC_SUCCESS);
    462: }

    464: /*@
    465:   DMDASetStencilType - Sets the type of the communication stencil

    467:   Logically Collective

    469:   Input Parameters:
    470: + da    - The `DMDA`
    471: - stype - The stencil type, use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.

    473:   Level: intermediate

    475: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMDAStencilType`, `DMDA_STENCIL_BOX`, `DMDA_STENCIL_STAR.`
    476: @*/
    477: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
    478: {
    479:   DM_DA *dd = (DM_DA *)da->data;

    481:   PetscFunctionBegin;
    484:   PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    485:   dd->stencil_type = stype;
    486:   PetscFunctionReturn(PETSC_SUCCESS);
    487: }

    489: /*@
    490:   DMDAGetStencilType - Gets the type of the communication stencil

    492:   Not Collective

    494:   Input Parameter:
    495: . da - The `DMDA`

    497:   Output Parameter:
    498: . stype - The stencil type, use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.

    500:   Level: intermediate

    502: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMDAStencilType`, `DMDA_STENCIL_BOX`, `DMDA_STENCIL_STAR.`
    503: @*/
    504: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
    505: {
    506:   DM_DA *dd = (DM_DA *)da->data;

    508:   PetscFunctionBegin;
    510:   PetscAssertPointer(stype, 2);
    511:   *stype = dd->stencil_type;
    512:   PetscFunctionReturn(PETSC_SUCCESS);
    513: }

    515: /*@
    516:   DMDASetStencilWidth - Sets the width of the communication stencil

    518:   Logically Collective

    520:   Input Parameters:
    521: + da    - The `DMDA`
    522: - width - The stencil width

    524:   Level: intermediate

    526: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMDAStencilType`, `DMDA_STENCIL_BOX`, `DMDA_STENCIL_STAR.`
    527: @*/
    528: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
    529: {
    530:   DM_DA *dd = (DM_DA *)da->data;

    532:   PetscFunctionBegin;
    535:   PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    536:   dd->s = width;
    537:   PetscFunctionReturn(PETSC_SUCCESS);
    538: }

    540: /*@
    541:   DMDAGetStencilWidth - Gets the width of the communication stencil

    543:   Not Collective

    545:   Input Parameter:
    546: . da - The `DMDA`

    548:   Output Parameter:
    549: . width - The stencil width

    551:   Level: intermediate

    553: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMDAStencilType`, `DMDA_STENCIL_BOX`, `DMDA_STENCIL_STAR.`
    554: @*/
    555: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
    556: {
    557:   DM_DA *dd = (DM_DA *)da->data;

    559:   PetscFunctionBegin;
    561:   PetscAssertPointer(width, 2);
    562:   *width = dd->s;
    563:   PetscFunctionReturn(PETSC_SUCCESS);
    564: }

    566: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da, PetscInt M, PetscInt m, const PetscInt lx[])
    567: {
    568:   PetscInt i, sum;

    570:   PetscFunctionBegin;
    571:   PetscCheck(M >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Global dimension not set");
    572:   for (i = sum = 0; i < m; i++) sum += lx[i];
    573:   PetscCheck(sum == M, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "Ownership ranges sum to %" PetscInt_FMT " but global dimension is %" PetscInt_FMT, sum, M);
    574:   PetscFunctionReturn(PETSC_SUCCESS);
    575: }

    577: /*@
    578:   DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process

    580:   Logically Collective

    582:   Input Parameters:
    583: + da - The `DMDA`
    584: . lx - array containing number of nodes in the X direction on each process, or `NULL`. If non-null, must be of length da->m
    585: . ly - array containing number of nodes in the Y direction on each process, or `NULL`. If non-null, must be of length da->n
    586: - lz - array containing number of nodes in the Z direction on each process, or `NULL`. If non-null, must be of length da->p.

    588:   Level: intermediate

    590:   Note:
    591:   These numbers are NOT multiplied by the number of dof per node.

    593: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`
    594: @*/
    595: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
    596: {
    597:   DM_DA *dd = (DM_DA *)da->data;

    599:   PetscFunctionBegin;
    601:   PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
    602:   if (lx) {
    603:     PetscCheck(dd->m >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Cannot set ownership ranges before setting number of processes");
    604:     PetscCall(DMDACheckOwnershipRanges_Private(da, dd->M, dd->m, lx));
    605:     if (!dd->lx) PetscCall(PetscMalloc1(dd->m, &dd->lx));
    606:     PetscCall(PetscArraycpy(dd->lx, lx, dd->m));
    607:   }
    608:   if (ly) {
    609:     PetscCheck(dd->n >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Cannot set ownership ranges before setting number of processes");
    610:     PetscCall(DMDACheckOwnershipRanges_Private(da, dd->N, dd->n, ly));
    611:     if (!dd->ly) PetscCall(PetscMalloc1(dd->n, &dd->ly));
    612:     PetscCall(PetscArraycpy(dd->ly, ly, dd->n));
    613:   }
    614:   if (lz) {
    615:     PetscCheck(dd->p >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Cannot set ownership ranges before setting number of processes");
    616:     PetscCall(DMDACheckOwnershipRanges_Private(da, dd->P, dd->p, lz));
    617:     if (!dd->lz) PetscCall(PetscMalloc1(dd->p, &dd->lz));
    618:     PetscCall(PetscArraycpy(dd->lz, lz, dd->p));
    619:   }
    620:   PetscFunctionReturn(PETSC_SUCCESS);
    621: }

    623: /*@
    624:   DMDASetInterpolationType - Sets the type of interpolation that will be
    625:   returned by `DMCreateInterpolation()`

    627:   Logically Collective

    629:   Input Parameters:
    630: + da    - initial distributed array
    631: - ctype - `DMDA_Q1` and `DMDA_Q0` are currently the only supported forms

    633:   Level: intermediate

    635:   Note:
    636:   You should call this on the coarser of the two `DMDA` you pass to `DMCreateInterpolation()`

    638: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMDestroy()`, `DMDAInterpolationType`,
    639:           `DMDA_Q1`, `DMDA_Q0`
    640: @*/
    641: PetscErrorCode DMDASetInterpolationType(DM da, DMDAInterpolationType ctype)
    642: {
    643:   DM_DA *dd = (DM_DA *)da->data;

    645:   PetscFunctionBegin;
    648:   dd->interptype = ctype;
    649:   PetscFunctionReturn(PETSC_SUCCESS);
    650: }

    652: /*@
    653:   DMDAGetInterpolationType - Gets the type of interpolation that will be
    654:   used by `DMCreateInterpolation()`

    656:   Not Collective

    658:   Input Parameter:
    659: . da - distributed array

    661:   Output Parameter:
    662: . ctype - interpolation type (`DMDA_Q1` and `DMDA_Q0` are currently the only supported forms)

    664:   Level: intermediate

    666: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAInterpolationType`, `DMDASetInterpolationType()`, `DMCreateInterpolation()`,
    667:           `DMDA_Q1`, `DMDA_Q0`
    668: @*/
    669: PetscErrorCode DMDAGetInterpolationType(DM da, DMDAInterpolationType *ctype)
    670: {
    671:   DM_DA *dd = (DM_DA *)da->data;

    673:   PetscFunctionBegin;
    675:   PetscAssertPointer(ctype, 2);
    676:   *ctype = dd->interptype;
    677:   PetscFunctionReturn(PETSC_SUCCESS);
    678: }

    680: /*@C
    681:   DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
    682:   processes neighbors.

    684:   Not Collective

    686:   Input Parameter:
    687: . da - the `DMDA` object

    689:   Output Parameter:
    690: . ranks - the neighbors ranks, stored with the x index increasing most rapidly. The process itself is in the list

    692:   Level: intermediate

    694:   Notes:
    695:   In 2d the `ranks` is of length 9, in 3d of length 27

    697:   Not supported in 1d

    699:   Do not free the array, it is freed when the `DMDA` is destroyed.

    701:   Fortran Note:
    702:   Use `DMDARestoreNeighbors()` to return the array when no longer needed

    704: .seealso: [](sec_struct), `DMDA`, `DM`
    705: @*/
    706: PetscErrorCode DMDAGetNeighbors(DM da, const PetscMPIInt *ranks[])
    707: {
    708:   DM_DA *dd = (DM_DA *)da->data;

    710:   PetscFunctionBegin;
    712:   *ranks = dd->neighbors;
    713:   PetscFunctionReturn(PETSC_SUCCESS);
    714: }

    716: /*@C
    717:   DMDAGetOwnershipRanges - Gets the number of indices in the x, y and z direction that are owned by each process in that direction

    719:   Not Collective

    721:   Input Parameter:
    722: . da - the `DMDA` object

    724:   Output Parameters:
    725: + lx - ownership along x direction (optional), its length is `m` the number of processes in the x-direction
    726: . ly - ownership along y direction (optional), its length is `n` the number of processes in the y-direction
    727: - lz - ownership along z direction (optional), its length is `p` the number of processes in the z-direction

    729:   Level: intermediate

    731:   Notes:
    732:   These correspond to the optional final arguments passed to `DMDACreate()`, `DMDACreate2d()`, `DMDACreate3d()`

    734:   You should not free these arrays, nor change the values in them. They will only have valid values while the
    735:   `DMDA` they came from still exists (has not been destroyed).

    737:   These numbers are NOT multiplied by the number of dof per node.

    739:   The meaning of these is different than that returned by `VecGetOwnerShipRanges()`

    741:   Fortran Notes:
    742:   Pass `PETSC_NULL_INT_POINTER` for any array not needed.

    744:   Use `DMDARestoreOwershipRange()` to return the arrays when no longer needed

    746: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetCorners()`, `DMDAGetGhostCorners()`, `DMDACreate()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `VecGetOwnershipRanges()`
    747: @*/
    748: PetscErrorCode DMDAGetOwnershipRanges(DM da, PeOp const PetscInt *lx[], PeOp const PetscInt *ly[], PeOp const PetscInt *lz[])
    749: {
    750:   DM_DA *dd = (DM_DA *)da->data;

    752:   PetscFunctionBegin;
    754:   if (lx) *lx = dd->lx;
    755:   if (ly) *ly = dd->ly;
    756:   if (lz) *lz = dd->lz;
    757:   PetscFunctionReturn(PETSC_SUCCESS);
    758: }

    760: /*@
    761:   DMDASetRefinementFactor - Set the ratios that the `DMDA` grid is refined

    763:   Logically Collective

    765:   Input Parameters:
    766: + da       - the `DMDA` object
    767: . refine_x - ratio of fine grid to coarse in the x direction (2 by default)
    768: . refine_y - ratio of fine grid to coarse in the y direction (2 by default)
    769: - refine_z - ratio of fine grid to coarse in the z direction (2 by default)

    771:   Options Database Keys:
    772: + -da_refine_x refine_x - refinement ratio in the x direction
    773: . -da_refine_y rafine_y - refinement ratio in the y direction
    774: . -da_refine_z refine_z - refinement ratio in the z direction
    775: - -da_refine n          - refine the `DMDA` object `n` times when it is created.

    777:   Level: intermediate

    779:   Note:
    780:   Pass `PETSC_IGNORE` to leave a value unchanged

    782: .seealso: [](sec_struct), `DM`, `DMDA`, `DMRefine()`, `DMDAGetRefinementFactor()`
    783: @*/
    784: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y, PetscInt refine_z)
    785: {
    786:   DM_DA *dd = (DM_DA *)da->data;

    788:   PetscFunctionBegin;

    794:   if (refine_x > 0) dd->refine_x = refine_x;
    795:   if (refine_y > 0) dd->refine_y = refine_y;
    796:   if (refine_z > 0) dd->refine_z = refine_z;
    797:   PetscFunctionReturn(PETSC_SUCCESS);
    798: }

    800: /*@
    801:   DMDAGetRefinementFactor - Gets the ratios that the `DMDA` grid is refined

    803:   Not Collective

    805:   Input Parameter:
    806: . da - the `DMDA` object

    808:   Output Parameters:
    809: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
    810: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
    811: - refine_z - ratio of fine grid to coarse in z direction (2 by default)

    813:   Level: intermediate

    815:   Note:
    816:   Pass `NULL` for values you do not need

    818: .seealso: [](sec_struct), `DM`, `DMDA`, `DMRefine()`, `DMDASetRefinementFactor()`
    819: @*/
    820: PetscErrorCode DMDAGetRefinementFactor(DM da, PeOp PetscInt *refine_x, PeOp PetscInt *refine_y, PeOp PetscInt *refine_z)
    821: {
    822:   DM_DA *dd = (DM_DA *)da->data;

    824:   PetscFunctionBegin;
    826:   if (refine_x) *refine_x = dd->refine_x;
    827:   if (refine_y) *refine_y = dd->refine_y;
    828:   if (refine_z) *refine_z = dd->refine_z;
    829:   PetscFunctionReturn(PETSC_SUCCESS);
    830: }

    832: /*@C
    833:   DMDASetGetMatrix - Sets the routine used by the `DMDA` to allocate a matrix.

    835:   Logically Collective; No Fortran Support

    837:   Input Parameters:
    838: + da - the `DMDA` object
    839: - f  - the function that allocates the matrix for that specific `DMDA`

    841:   Calling sequence of `f`:
    842: + da - the `DMDA` object
    843: - A  - the created matrix

    845:   Level: developer

    847:   Notes:
    848:   If the function is not provided a default function is used that uses the `DMDAStencilType`, `DMBoundaryType`, and value of `DMDASetStencilWidth()`
    849:   to construct the matrix.

    851:   See `DMDASetBlockFills()` that provides a simple way to provide the nonzero structure for
    852:   the diagonal and off-diagonal blocks of the matrix without providing a custom function

    854:   Developer Note:
    855:   This should be called `DMDASetCreateMatrix()`

    857: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateMatrix()`, `DMDASetBlockFills()`
    858: @*/
    859: PetscErrorCode DMDASetGetMatrix(DM da, PetscErrorCode (*f)(DM da, Mat *A))
    860: {
    861:   PetscFunctionBegin;
    863:   da->ops->creatematrix = f;
    864:   PetscFunctionReturn(PETSC_SUCCESS);
    865: }

    867: /*@
    868:   DMDAMapMatStencilToGlobal - Map a list of `MatStencil` on a grid to global indices.

    870:   Not Collective

    872:   Input Parameters:
    873: + da   - the `DMDA` object
    874: . m    - number of `MatStencil` to map
    875: - idxm - grid points (and component number when dof > 1)

    877:   Output Parameter:
    878: . gidxm - global row indices

    880:   Level: intermediate

    882: .seealso: [](sec_struct), `DM`, `DMDA`, `MatStencil`
    883: @*/
    884: PetscErrorCode DMDAMapMatStencilToGlobal(DM da, PetscInt m, const MatStencil idxm[], PetscInt gidxm[])
    885: {
    886:   const DM_DA           *dd  = (const DM_DA *)da->data;
    887:   const PetscInt        *dxm = (const PetscInt *)idxm;
    888:   PetscInt               i, j, sdim, tmp, dim;
    889:   PetscInt               dims[4], starts[4], dims2[3], starts2[3], dof = dd->w;
    890:   ISLocalToGlobalMapping ltog;

    892:   PetscFunctionBegin;
    893:   if (m <= 0) PetscFunctionReturn(PETSC_SUCCESS);

    895:   /* Code adapted from DMDAGetGhostCorners() */
    896:   starts2[0] = dd->Xs / dof + dd->xo;
    897:   starts2[1] = dd->Ys + dd->yo;
    898:   starts2[2] = dd->Zs + dd->zo;
    899:   dims2[0]   = (dd->Xe - dd->Xs) / dof;
    900:   dims2[1]   = (dd->Ye - dd->Ys);
    901:   dims2[2]   = (dd->Ze - dd->Zs);

    903:   /* As if we do MatSetStencil() to get dims[]/starts[] of mat->stencil */
    904:   dim  = da->dim;                   /* DA dim: 1 to 3 */
    905:   sdim = dim + (dof > 1 ? 1 : 0);   /* Dimensions in MatStencil's (k,j,i,c) view */
    906:   for (i = 0; i < dim; i++) {       /* Reverse the order and also skip the unused dimensions */
    907:     dims[i]   = dims2[dim - i - 1]; /* ex. dims/starts[] are in order of {i} for 1D, {j,i} for 2D and {k,j,i} for 3D */
    908:     starts[i] = starts2[dim - i - 1];
    909:   }
    910:   starts[dim] = 0; /* Append the extra dim for dof (won't be used below if dof=1) */
    911:   dims[dim]   = dof;

    913:   /* Map stencils to local indices (code adapted from MatSetValuesStencil()) */
    914:   for (i = 0; i < m; i++) {
    915:     dxm += 3 - dim; /* Input is {k,j,i,c}; move the pointer to the first used index, e.g., j in 2D */
    916:     tmp = 0;
    917:     for (j = 0; j < sdim; j++) {                                                      /* Iter over, ex. j,i or j,i,c in 2D */
    918:       if (tmp < 0 || dxm[j] < starts[j] || dxm[j] >= (starts[j] + dims[j])) tmp = -1; /* Beyond the ghost region, therefore ignored with negative indices */
    919:       else tmp = tmp * dims[j] + (dxm[j] - starts[j]);
    920:     }
    921:     gidxm[i] = tmp;
    922:     /* Move to the next MatStencil point */
    923:     if (dof > 1) dxm += sdim; /* c is already counted in sdim */
    924:     else dxm += sdim + 1;     /* skip the unused c */
    925:   }

    927:   /* Map local indices to global indices */
    928:   PetscCall(DMGetLocalToGlobalMapping(da, &ltog));
    929:   PetscCall(ISLocalToGlobalMappingApply(ltog, m, gidxm, gidxm));
    930:   PetscFunctionReturn(PETSC_SUCCESS);
    931: }

    933: /*
    934:   Creates "balanced" ownership ranges after refinement, constrained by the need for the
    935:   fine grid boundaries to fall within one stencil width of the coarse partition.

    937:   Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
    938: */
    939: static PetscErrorCode DMDARefineOwnershipRanges(DM da, PetscBool periodic, PetscInt stencil_width, PetscInt ratio, PetscInt m, const PetscInt lc[], PetscInt lf[])
    940: {
    941:   PetscInt i, totalc = 0, remaining, startc = 0, startf = 0;

    943:   PetscFunctionBegin;
    944:   PetscCheck(ratio >= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Requested refinement ratio %" PetscInt_FMT " must be at least 1", ratio);
    945:   if (ratio == 1) {
    946:     PetscCall(PetscArraycpy(lf, lc, m));
    947:     PetscFunctionReturn(PETSC_SUCCESS);
    948:   }
    949:   for (i = 0; i < m; i++) totalc += lc[i];
    950:   remaining = (!periodic) + ratio * (totalc - (!periodic));
    951:   for (i = 0; i < m; i++) {
    952:     PetscInt want = remaining / (m - i) + !!(remaining % (m - i));
    953:     if (i == m - 1) lf[i] = want;
    954:     else {
    955:       const PetscInt nextc = startc + lc[i];
    956:       /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
    957:        * coarse stencil width of the first coarse node in the next subdomain. */
    958:       while ((startf + want) / ratio < nextc - stencil_width) want++;
    959:       /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
    960:        * coarse stencil width of the last coarse node in the current subdomain. */
    961:       while ((startf + want - 1 + ratio - 1) / ratio > nextc - 1 + stencil_width) want--;
    962:       /* Make sure all constraints are satisfied */
    963:       if (want < 0 || want > remaining || ((startf + want) / ratio < nextc - stencil_width) || ((startf + want - 1 + ratio - 1) / ratio > nextc - 1 + stencil_width))
    964:         SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Could not find a compatible refined ownership range");
    965:     }
    966:     lf[i] = want;
    967:     startc += lc[i];
    968:     startf += lf[i];
    969:     remaining -= lf[i];
    970:   }
    971:   PetscFunctionReturn(PETSC_SUCCESS);
    972: }

    974: /*
    975:   Creates "balanced" ownership ranges after coarsening, constrained by the need for the
    976:   fine grid boundaries to fall within one stencil width of the coarse partition.

    978:   Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
    979: */
    980: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da, PetscBool periodic, PetscInt stencil_width, PetscInt ratio, PetscInt m, const PetscInt lf[], PetscInt lc[])
    981: {
    982:   PetscInt i, totalf, remaining, startc, startf;

    984:   PetscFunctionBegin;
    985:   PetscCheck(ratio >= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Requested refinement ratio %" PetscInt_FMT " must be at least 1", ratio);
    986:   if (ratio == 1) {
    987:     PetscCall(PetscArraycpy(lc, lf, m));
    988:     PetscFunctionReturn(PETSC_SUCCESS);
    989:   }
    990:   for (i = 0, totalf = 0; i < m; i++) totalf += lf[i];
    991:   remaining = (!periodic) + (totalf - (!periodic)) / ratio;
    992:   for (i = 0, startc = 0, startf = 0; i < m; i++) {
    993:     PetscInt want = remaining / (m - i) + !!(remaining % (m - i));
    994:     if (i == m - 1) lc[i] = want;
    995:     else {
    996:       const PetscInt nextf = startf + lf[i];
    997:       /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
    998:        * node is within one stencil width. */
    999:       while (nextf / ratio < startc + want - stencil_width) want--;
    1000:       /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
    1001:        * fine node is within one stencil width. */
    1002:       while ((nextf - 1 + ratio - 1) / ratio > startc + want - 1 + stencil_width) want++;
    1003:       if (want < 0 || want > remaining || (nextf / ratio < startc + want - stencil_width) || ((nextf - 1 + ratio - 1) / ratio > startc + want - 1 + stencil_width))
    1004:         SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Could not find a compatible coarsened ownership range");
    1005:     }
    1006:     lc[i] = want;
    1007:     startc += lc[i];
    1008:     startf += lf[i];
    1009:     remaining -= lc[i];
    1010:   }
    1011:   PetscFunctionReturn(PETSC_SUCCESS);
    1012: }

    1014: PetscErrorCode DMRefine_DA(DM da, MPI_Comm comm, DM *daref)
    1015: {
    1016:   PetscInt M, N, P, i, dim;
    1017:   Vec      coordsc, coordsf;
    1018:   DM       da2;
    1019:   DM_DA   *dd = (DM_DA *)da->data, *dd2;

    1021:   PetscFunctionBegin;
    1023:   PetscAssertPointer(daref, 3);

    1025:   PetscCall(DMGetDimension(da, &dim));
    1026:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
    1027:     M = dd->refine_x * dd->M;
    1028:   } else {
    1029:     M = 1 + dd->refine_x * (dd->M - 1);
    1030:   }
    1031:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
    1032:     if (dim > 1) {
    1033:       N = dd->refine_y * dd->N;
    1034:     } else {
    1035:       N = 1;
    1036:     }
    1037:   } else {
    1038:     N = 1 + dd->refine_y * (dd->N - 1);
    1039:   }
    1040:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
    1041:     if (dim > 2) {
    1042:       P = dd->refine_z * dd->P;
    1043:     } else {
    1044:       P = 1;
    1045:     }
    1046:   } else {
    1047:     P = 1 + dd->refine_z * (dd->P - 1);
    1048:   }
    1049:   PetscCall(DMDACreate(PetscObjectComm((PetscObject)da), &da2));
    1050:   PetscCall(DMSetOptionsPrefix(da2, ((PetscObject)da)->prefix));
    1051:   PetscCall(DMSetDimension(da2, dim));
    1052:   PetscCall(DMDASetSizes(da2, M, N, P));
    1053:   PetscCall(DMDASetNumProcs(da2, dd->m, dd->n, dd->p));
    1054:   PetscCall(DMDASetBoundaryType(da2, dd->bx, dd->by, dd->bz));
    1055:   PetscCall(DMDASetDof(da2, dd->w));
    1056:   PetscCall(DMDASetStencilType(da2, dd->stencil_type));
    1057:   PetscCall(DMDASetStencilWidth(da2, dd->s));
    1058:   if (dim == 3) {
    1059:     PetscInt *lx, *ly, *lz;
    1060:     PetscCall(PetscMalloc3(dd->m, &lx, dd->n, &ly, dd->p, &lz));
    1061:     PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx));
    1062:     PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_y, dd->n, dd->ly, ly));
    1063:     PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_z, dd->p, dd->lz, lz));
    1064:     PetscCall(DMDASetOwnershipRanges(da2, lx, ly, lz));
    1065:     PetscCall(PetscFree3(lx, ly, lz));
    1066:   } else if (dim == 2) {
    1067:     PetscInt *lx, *ly;
    1068:     PetscCall(PetscMalloc2(dd->m, &lx, dd->n, &ly));
    1069:     PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx));
    1070:     PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_y, dd->n, dd->ly, ly));
    1071:     PetscCall(DMDASetOwnershipRanges(da2, lx, ly, NULL));
    1072:     PetscCall(PetscFree2(lx, ly));
    1073:   } else if (dim == 1) {
    1074:     PetscInt *lx;
    1075:     PetscCall(PetscMalloc1(dd->m, &lx));
    1076:     PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx));
    1077:     PetscCall(DMDASetOwnershipRanges(da2, lx, NULL, NULL));
    1078:     PetscCall(PetscFree(lx));
    1079:   }
    1080:   dd2 = (DM_DA *)da2->data;

    1082:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
    1083:   da2->ops->creatematrix = da->ops->creatematrix;
    1084:   /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
    1085:   da2->ops->getcoloring = da->ops->getcoloring;
    1086:   dd2->interptype       = dd->interptype;

    1088:   /* copy fill information if given */
    1089:   if (dd->dfill) {
    1090:     PetscCall(PetscMalloc1(dd->dfill[dd->w] + dd->w + 1, &dd2->dfill));
    1091:     PetscCall(PetscArraycpy(dd2->dfill, dd->dfill, dd->dfill[dd->w] + dd->w + 1));
    1092:   }
    1093:   if (dd->ofill) {
    1094:     PetscCall(PetscMalloc1(dd->ofill[dd->w] + dd->w + 1, &dd2->ofill));
    1095:     PetscCall(PetscArraycpy(dd2->ofill, dd->ofill, dd->ofill[dd->w] + dd->w + 1));
    1096:   }
    1097:   /* copy the refine information */
    1098:   dd2->coarsen_x = dd2->refine_x = dd->refine_x;
    1099:   dd2->coarsen_y = dd2->refine_y = dd->refine_y;
    1100:   dd2->coarsen_z = dd2->refine_z = dd->refine_z;

    1102:   if (dd->refine_z_hier) {
    1103:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
    1104:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
    1105:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
    1106:     PetscCall(PetscMalloc1(dd2->refine_z_hier_n, &dd2->refine_z_hier));
    1107:     PetscCall(PetscArraycpy(dd2->refine_z_hier, dd->refine_z_hier, dd2->refine_z_hier_n));
    1108:   }
    1109:   if (dd->refine_y_hier) {
    1110:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
    1111:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
    1112:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
    1113:     PetscCall(PetscMalloc1(dd2->refine_y_hier_n, &dd2->refine_y_hier));
    1114:     PetscCall(PetscArraycpy(dd2->refine_y_hier, dd->refine_y_hier, dd2->refine_y_hier_n));
    1115:   }
    1116:   if (dd->refine_x_hier) {
    1117:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
    1118:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
    1119:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
    1120:     PetscCall(PetscMalloc1(dd2->refine_x_hier_n, &dd2->refine_x_hier));
    1121:     PetscCall(PetscArraycpy(dd2->refine_x_hier, dd->refine_x_hier, dd2->refine_x_hier_n));
    1122:   }

    1124:   /* copy vector type information */
    1125:   PetscCall(DMSetVecType(da2, da->vectype));

    1127:   dd2->lf = dd->lf;
    1128:   dd2->lj = dd->lj;

    1130:   da2->leveldown = da->leveldown;
    1131:   da2->levelup   = da->levelup + 1;

    1133:   PetscCall(DMSetUp(da2));

    1135:   /* interpolate coordinates if they are set on the coarse grid */
    1136:   PetscCall(DMGetCoordinates(da, &coordsc));
    1137:   if (coordsc) {
    1138:     DM  cdaf, cdac;
    1139:     Mat II;

    1141:     PetscCall(DMGetCoordinateDM(da, &cdac));
    1142:     PetscCall(DMGetCoordinateDM(da2, &cdaf));
    1143:     /* force creation of the coordinate vector */
    1144:     PetscCall(DMDASetUniformCoordinates(da2, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
    1145:     PetscCall(DMGetCoordinates(da2, &coordsf));
    1146:     PetscCall(DMCreateInterpolation(cdac, cdaf, &II, NULL));
    1147:     PetscCall(MatInterpolate(II, coordsc, coordsf));
    1148:     PetscCall(MatDestroy(&II));
    1149:   }

    1151:   for (i = 0; i < da->bs; i++) {
    1152:     const char *fieldname;
    1153:     PetscCall(DMDAGetFieldName(da, i, &fieldname));
    1154:     PetscCall(DMDASetFieldName(da2, i, fieldname));
    1155:   }

    1157:   *daref = da2;
    1158:   PetscFunctionReturn(PETSC_SUCCESS);
    1159: }

    1161: PetscErrorCode DMCoarsen_DA(DM dmf, MPI_Comm comm, DM *dmc)
    1162: {
    1163:   PetscInt M, N, P, i, dim;
    1164:   Vec      coordsc, coordsf;
    1165:   DM       dmc2;
    1166:   DM_DA   *dd = (DM_DA *)dmf->data, *dd2;

    1168:   PetscFunctionBegin;
    1170:   PetscAssertPointer(dmc, 3);

    1172:   PetscCall(DMGetDimension(dmf, &dim));
    1173:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
    1174:     M = dd->M / dd->coarsen_x;
    1175:   } else {
    1176:     M = 1 + (dd->M - 1) / dd->coarsen_x;
    1177:   }
    1178:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
    1179:     if (dim > 1) {
    1180:       N = dd->N / dd->coarsen_y;
    1181:     } else {
    1182:       N = 1;
    1183:     }
    1184:   } else {
    1185:     N = 1 + (dd->N - 1) / dd->coarsen_y;
    1186:   }
    1187:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
    1188:     if (dim > 2) {
    1189:       P = dd->P / dd->coarsen_z;
    1190:     } else {
    1191:       P = 1;
    1192:     }
    1193:   } else {
    1194:     P = 1 + (dd->P - 1) / dd->coarsen_z;
    1195:   }
    1196:   PetscCall(DMDACreate(PetscObjectComm((PetscObject)dmf), &dmc2));
    1197:   PetscCall(DMSetOptionsPrefix(dmc2, ((PetscObject)dmf)->prefix));
    1198:   PetscCall(DMSetDimension(dmc2, dim));
    1199:   PetscCall(DMDASetSizes(dmc2, M, N, P));
    1200:   PetscCall(DMDASetNumProcs(dmc2, dd->m, dd->n, dd->p));
    1201:   PetscCall(DMDASetBoundaryType(dmc2, dd->bx, dd->by, dd->bz));
    1202:   PetscCall(DMDASetDof(dmc2, dd->w));
    1203:   PetscCall(DMDASetStencilType(dmc2, dd->stencil_type));
    1204:   PetscCall(DMDASetStencilWidth(dmc2, dd->s));
    1205:   if (dim == 3) {
    1206:     PetscInt *lx, *ly, *lz;
    1207:     PetscCall(PetscMalloc3(dd->m, &lx, dd->n, &ly, dd->p, &lz));
    1208:     PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx));
    1209:     PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_y, dd->n, dd->ly, ly));
    1210:     PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_z, dd->p, dd->lz, lz));
    1211:     PetscCall(DMDASetOwnershipRanges(dmc2, lx, ly, lz));
    1212:     PetscCall(PetscFree3(lx, ly, lz));
    1213:   } else if (dim == 2) {
    1214:     PetscInt *lx, *ly;
    1215:     PetscCall(PetscMalloc2(dd->m, &lx, dd->n, &ly));
    1216:     PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx));
    1217:     PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_y, dd->n, dd->ly, ly));
    1218:     PetscCall(DMDASetOwnershipRanges(dmc2, lx, ly, NULL));
    1219:     PetscCall(PetscFree2(lx, ly));
    1220:   } else if (dim == 1) {
    1221:     PetscInt *lx;
    1222:     PetscCall(PetscMalloc1(dd->m, &lx));
    1223:     PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx));
    1224:     PetscCall(DMDASetOwnershipRanges(dmc2, lx, NULL, NULL));
    1225:     PetscCall(PetscFree(lx));
    1226:   }
    1227:   dd2 = (DM_DA *)dmc2->data;

    1229:   /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
    1230:   /* dmc2->ops->createinterpolation = dmf->ops->createinterpolation; copying this one causes trouble for DMSetVI */
    1231:   dmc2->ops->creatematrix = dmf->ops->creatematrix;
    1232:   dmc2->ops->getcoloring  = dmf->ops->getcoloring;
    1233:   dd2->interptype         = dd->interptype;

    1235:   /* copy fill information if given */
    1236:   if (dd->dfill) {
    1237:     PetscCall(PetscMalloc1(dd->dfill[dd->w] + dd->w + 1, &dd2->dfill));
    1238:     PetscCall(PetscArraycpy(dd2->dfill, dd->dfill, dd->dfill[dd->w] + dd->w + 1));
    1239:   }
    1240:   if (dd->ofill) {
    1241:     PetscCall(PetscMalloc1(dd->ofill[dd->w] + dd->w + 1, &dd2->ofill));
    1242:     PetscCall(PetscArraycpy(dd2->ofill, dd->ofill, dd->ofill[dd->w] + dd->w + 1));
    1243:   }
    1244:   /* copy the refine information */
    1245:   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
    1246:   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
    1247:   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;

    1249:   if (dd->refine_z_hier) {
    1250:     if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_z_hier_n) dd2->refine_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 1];
    1251:     if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_z_hier_n) dd2->coarsen_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 2];
    1252:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
    1253:     PetscCall(PetscMalloc1(dd2->refine_z_hier_n, &dd2->refine_z_hier));
    1254:     PetscCall(PetscArraycpy(dd2->refine_z_hier, dd->refine_z_hier, dd2->refine_z_hier_n));
    1255:   }
    1256:   if (dd->refine_y_hier) {
    1257:     if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_y_hier_n) dd2->refine_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 1];
    1258:     if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_y_hier_n) dd2->coarsen_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 2];
    1259:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
    1260:     PetscCall(PetscMalloc1(dd2->refine_y_hier_n, &dd2->refine_y_hier));
    1261:     PetscCall(PetscArraycpy(dd2->refine_y_hier, dd->refine_y_hier, dd2->refine_y_hier_n));
    1262:   }
    1263:   if (dd->refine_x_hier) {
    1264:     if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_x_hier_n) dd2->refine_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 1];
    1265:     if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_x_hier_n) dd2->coarsen_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 2];
    1266:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
    1267:     PetscCall(PetscMalloc1(dd2->refine_x_hier_n, &dd2->refine_x_hier));
    1268:     PetscCall(PetscArraycpy(dd2->refine_x_hier, dd->refine_x_hier, dd2->refine_x_hier_n));
    1269:   }

    1271:   /* copy vector type information */
    1272:   PetscCall(DMSetVecType(dmc2, dmf->vectype));

    1274:   dd2->lf = dd->lf;
    1275:   dd2->lj = dd->lj;

    1277:   dmc2->leveldown = dmf->leveldown + 1;
    1278:   dmc2->levelup   = dmf->levelup;

    1280:   PetscCall(DMSetUp(dmc2));

    1282:   /* inject coordinates if they are set on the fine grid */
    1283:   PetscCall(DMGetCoordinates(dmf, &coordsf));
    1284:   if (coordsf) {
    1285:     DM         cdaf, cdac;
    1286:     Mat        inject;
    1287:     VecScatter vscat;

    1289:     PetscCall(DMGetCoordinateDM(dmf, &cdaf));
    1290:     PetscCall(DMGetCoordinateDM(dmc2, &cdac));
    1291:     /* force creation of the coordinate vector */
    1292:     PetscCall(DMDASetUniformCoordinates(dmc2, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
    1293:     PetscCall(DMGetCoordinates(dmc2, &coordsc));

    1295:     PetscCall(DMCreateInjection(cdac, cdaf, &inject));
    1296:     PetscCall(MatScatterGetVecScatter(inject, &vscat));
    1297:     PetscCall(VecScatterBegin(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD));
    1298:     PetscCall(VecScatterEnd(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD));
    1299:     PetscCall(MatDestroy(&inject));
    1300:   }

    1302:   for (i = 0; i < dmf->bs; i++) {
    1303:     const char *fieldname;
    1304:     PetscCall(DMDAGetFieldName(dmf, i, &fieldname));
    1305:     PetscCall(DMDASetFieldName(dmc2, i, fieldname));
    1306:   }

    1308:   *dmc = dmc2;
    1309:   PetscFunctionReturn(PETSC_SUCCESS);
    1310: }

    1312: PetscErrorCode DMRefineHierarchy_DA(DM da, PetscInt nlevels, DM daf[])
    1313: {
    1314:   PetscInt i, n, *refx, *refy, *refz;

    1316:   PetscFunctionBegin;
    1318:   PetscCheck(nlevels >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "nlevels cannot be negative");
    1319:   if (nlevels == 0) PetscFunctionReturn(PETSC_SUCCESS);
    1320:   PetscAssertPointer(daf, 3);

    1322:   /* Get refinement factors, defaults taken from the coarse DMDA */
    1323:   PetscCall(PetscMalloc3(nlevels, &refx, nlevels, &refy, nlevels, &refz));
    1324:   for (i = 0; i < nlevels; i++) PetscCall(DMDAGetRefinementFactor(da, &refx[i], &refy[i], &refz[i]));
    1325:   n = nlevels;
    1326:   PetscCall(PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_x", refx, &n, NULL));
    1327:   n = nlevels;
    1328:   PetscCall(PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_y", refy, &n, NULL));
    1329:   n = nlevels;
    1330:   PetscCall(PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_z", refz, &n, NULL));

    1332:   PetscCall(DMDASetRefinementFactor(da, refx[0], refy[0], refz[0]));
    1333:   PetscCall(DMRefine(da, PetscObjectComm((PetscObject)da), &daf[0]));
    1334:   for (i = 1; i < nlevels; i++) {
    1335:     PetscCall(DMDASetRefinementFactor(daf[i - 1], refx[i], refy[i], refz[i]));
    1336:     PetscCall(DMRefine(daf[i - 1], PetscObjectComm((PetscObject)da), &daf[i]));
    1337:   }
    1338:   PetscCall(PetscFree3(refx, refy, refz));
    1339:   PetscFunctionReturn(PETSC_SUCCESS);
    1340: }

    1342: PetscErrorCode DMCoarsenHierarchy_DA(DM da, PetscInt nlevels, DM dac[])
    1343: {
    1344:   PetscInt i;

    1346:   PetscFunctionBegin;
    1348:   PetscCheck(nlevels >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "nlevels cannot be negative");
    1349:   if (nlevels == 0) PetscFunctionReturn(PETSC_SUCCESS);
    1350:   PetscAssertPointer(dac, 3);
    1351:   PetscCall(DMCoarsen(da, PetscObjectComm((PetscObject)da), &dac[0]));
    1352:   for (i = 1; i < nlevels; i++) PetscCall(DMCoarsen(dac[i - 1], PetscObjectComm((PetscObject)da), &dac[i]));
    1353:   PetscFunctionReturn(PETSC_SUCCESS);
    1354: }

    1356: static PetscErrorCode DMDASetGLLCoordinates_1d(DM dm, PetscInt n, PetscReal *nodes)
    1357: {
    1358:   PetscInt     i, j, xs, xn, q;
    1359:   PetscScalar *xx;
    1360:   PetscReal    h;
    1361:   Vec          x;
    1362:   DM_DA       *da = (DM_DA *)dm->data;

    1364:   PetscFunctionBegin;
    1365:   PetscCheck(da->bx != DM_BOUNDARY_PERIODIC, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not yet implemented for periodic");
    1366:   PetscCall(DMDAGetInfo(dm, NULL, &q, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    1367:   q = (q - 1) / (n - 1); /* number of spectral elements */
    1368:   h = 2.0 / q;
    1369:   PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xn, NULL, NULL));
    1370:   xs = xs / (n - 1);
    1371:   xn = xn / (n - 1);
    1372:   PetscCall(DMDASetUniformCoordinates(dm, -1., 1., 0., 0., 0., 0.));
    1373:   PetscCall(DMGetCoordinates(dm, &x));
    1374:   PetscCall(DMDAVecGetArray(dm, x, &xx));

    1376:   /* loop over local spectral elements */
    1377:   for (j = xs; j < xs + xn; j++) {
    1378:     /*
    1379:      Except for the first process, each process starts on the second GLL point of the first element on that process
    1380:      */
    1381:     for (i = (j == xs && xs > 0) ? 1 : 0; i < n; i++) xx[j * (n - 1) + i] = -1.0 + h * j + h * (nodes[i] + 1.0) / 2.;
    1382:   }
    1383:   PetscCall(DMDAVecRestoreArray(dm, x, &xx));
    1384:   PetscFunctionReturn(PETSC_SUCCESS);
    1385: }

    1387: /*@
    1388:   DMDASetGLLCoordinates - Sets the global coordinates from -1 to 1 to the GLL points of as many GLL elements that fit the number of grid points

    1390:   Collective

    1392:   Input Parameters:
    1393: + da    - the `DMDA` object
    1394: . n     - the number of GLL nodes
    1395: - nodes - the GLL nodes

    1397:   Level: advanced

    1399:   Note:
    1400:   The parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
    1401:   on each process much be divisible by the number of GLL elements needed per process. This depends on whether the `DMDA` is
    1402:   periodic or not.

    1404: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `PetscDTGaussLobattoLegendreQuadrature()`, `DMGetCoordinates()`
    1405: @*/
    1406: PetscErrorCode DMDASetGLLCoordinates(DM da, PetscInt n, PetscReal *nodes)
    1407: {
    1408:   PetscFunctionBegin;
    1409:   PetscCheck(da->dim == 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not yet implemented for 2 or 3d");
    1410:   PetscCall(DMDASetGLLCoordinates_1d(da, n, nodes));
    1411:   PetscFunctionReturn(PETSC_SUCCESS);
    1412: }

    1414: PetscErrorCode DMGetCompatibility_DA(DM da1, DM dm2, PetscBool *compatible, PetscBool *set)
    1415: {
    1416:   DM_DA    *dd1 = (DM_DA *)da1->data, *dd2;
    1417:   DM        da2;
    1418:   DMType    dmtype2;
    1419:   PetscBool isda, compatibleLocal;
    1420:   PetscInt  i;

    1422:   PetscFunctionBegin;
    1423:   PetscCheck(da1->setupcalled, PetscObjectComm((PetscObject)da1), PETSC_ERR_ARG_WRONGSTATE, "DMSetUp() must be called on first DM before DMGetCompatibility()");
    1424:   PetscCall(DMGetType(dm2, &dmtype2));
    1425:   PetscCall(PetscStrcmp(dmtype2, DMDA, &isda));
    1426:   if (isda) {
    1427:     da2 = dm2;
    1428:     dd2 = (DM_DA *)da2->data;
    1429:     PetscCheck(da2->setupcalled, PetscObjectComm((PetscObject)da2), PETSC_ERR_ARG_WRONGSTATE, "DMSetUp() must be called on second DM before DMGetCompatibility()");
    1430:     compatibleLocal = (PetscBool)(da1->dim == da2->dim);
    1431:     if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
    1432:     /*                                                                           Global size              ranks               Boundary type */
    1433:     if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
    1434:     if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
    1435:     if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
    1436:     if (compatibleLocal) {
    1437:       for (i = 0; i < dd1->m; ++i) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i])); /* Local size     */
    1438:     }
    1439:     if (compatibleLocal && da1->dim > 1) {
    1440:       for (i = 0; i < dd1->n; ++i) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
    1441:     }
    1442:     if (compatibleLocal && da1->dim > 2) {
    1443:       for (i = 0; i < dd1->p; ++i) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
    1444:     }
    1445:     *compatible = compatibleLocal;
    1446:     *set        = PETSC_TRUE;
    1447:   } else {
    1448:     /* Decline to determine compatibility with other DM types */
    1449:     *set = PETSC_FALSE;
    1450:   }
    1451:   PetscFunctionReturn(PETSC_SUCCESS);
    1452: }