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 on da

  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: `DM`, `DMDA`, `PetscSplitOwnership()`
 20: @*/
 21: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
 22: {
 23:   DM_DA *dd = (DM_DA *)da->data;


 34:   dd->M = M;
 35:   dd->N = N;
 36:   dd->P = P;
 37:   return 0;
 38: }

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

 43:   Logically Collective on da

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

 51:   Level: intermediate

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

 64:   dd->m = m;
 65:   dd->n = n;
 66:   dd->p = p;
 67:   if (da->dim == 2) {
 68:     PetscMPIInt size;
 69:     MPI_Comm_size(PetscObjectComm((PetscObject)da), &size);
 70:     if ((dd->m > 0) && (dd->n < 0)) {
 71:       dd->n = size / dd->m;
 73:     }
 74:     if ((dd->n > 0) && (dd->m < 0)) {
 75:       dd->m = size / dd->n;
 77:     }
 78:   }
 79:   return 0;
 80: }

 82: /*@
 83:   DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.

 85:   Not collective

 87:   Input Parameters:
 88: + da    - The `DMDA`
 89: - bx,by,bz - One of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`

 91:   Level: intermediate

 93: .seealso: `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMBoundaryType`
 94: @*/
 95: PetscErrorCode DMDASetBoundaryType(DM da, DMBoundaryType bx, DMBoundaryType by, DMBoundaryType bz)
 96: {
 97:   DM_DA *dd = (DM_DA *)da->data;

104:   dd->bx = bx;
105:   dd->by = by;
106:   dd->bz = bz;
107:   return 0;
108: }

110: /*@
111:   DMDASetDof - Sets the number of degrees of freedom per vertex

113:   Not collective

115:   Input Parameters:
116: + da  - The `DMDA`
117: - dof - Number of degrees of freedom

119:   Level: intermediate

121: .seealso: `DM`, `DMDA`, `DMDAGetDof()`, `DMDACreate()`, `DMDestroy()`
122: @*/
123: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
124: {
125:   DM_DA *dd = (DM_DA *)da->data;

130:   dd->w  = dof;
131:   da->bs = dof;
132:   return 0;
133: }

135: /*@
136:   DMDAGetDof - Gets the number of degrees of freedom per vertex

138:   Not collective

140:   Input Parameter:
141: . da  - The `DMDA`

143:   Output Parameter:
144: . dof - Number of degrees of freedom

146:   Level: intermediate

148: .seealso: `DM`, `DMDA`, `DMDASetDof()`, `DMDACreate()`, `DMDestroy()`
149: @*/
150: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
151: {
152:   DM_DA *dd = (DM_DA *)da->data;

156:   *dof = dd->w;
157:   return 0;
158: }

160: /*@
161:   DMDAGetOverlap - Gets the size of the per-processor overlap.

163:   Not collective

165:   Input Parameter:
166: . da  - The `DMDA`

168:   Output Parameters:
169: + x   - Overlap in the x direction
170: . y   - Overlap in the y direction
171: - z   - Overlap in the z direction

173:   Level: intermediate

175: .seealso: `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDASetOverlap()`
176: @*/
177: PetscErrorCode DMDAGetOverlap(DM da, PetscInt *x, PetscInt *y, PetscInt *z)
178: {
179:   DM_DA *dd = (DM_DA *)da->data;

182:   if (x) *x = dd->xol;
183:   if (y) *y = dd->yol;
184:   if (z) *z = dd->zol;
185:   return 0;
186: }

188: /*@
189:   DMDASetOverlap - Sets the size of the per-processor overlap.

191:   Not collective

193:   Input Parameters:
194: + da  - The `DMDA`
195: . x   - Overlap in the x direction
196: . y   - Overlap in the y direction
197: - z   - Overlap in the z direction

199:   Level: intermediate

201: .seealso: `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDAGetOverlap()`
202: @*/
203: PetscErrorCode DMDASetOverlap(DM da, PetscInt x, PetscInt y, PetscInt z)
204: {
205:   DM_DA *dd = (DM_DA *)da->data;

211:   dd->xol = x;
212:   dd->yol = y;
213:   dd->zol = z;
214:   return 0;
215: }

217: /*@
218:   DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.

220:   Not collective

222:   Input Parameters:
223: . da  - The `DMDA`

225:   Output Parameters:
226: . Nsub   - Number of local subdomains created upon decomposition

228:   Level: intermediate

230: .seealso: `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDASetNumLocalSubDomains()`
231: @*/
232: PetscErrorCode DMDAGetNumLocalSubDomains(DM da, PetscInt *Nsub)
233: {
234:   DM_DA *dd = (DM_DA *)da->data;

237:   if (Nsub) *Nsub = dd->Nsub;
238:   return 0;
239: }

241: /*@
242:   DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.

244:   Not collective

246:   Input Parameters:
247: + da  - The `DMDA`
248: - Nsub - The number of local subdomains requested

250:   Level: intermediate

252: .seealso: `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDAGetNumLocalSubDomains()`
253: @*/
254: PetscErrorCode DMDASetNumLocalSubDomains(DM da, PetscInt Nsub)
255: {
256:   DM_DA *dd = (DM_DA *)da->data;

260:   dd->Nsub = Nsub;
261:   return 0;
262: }

264: /*@
265:   DMDASetOffset - Sets the index offset of the DA.

267:   Collective on da

269:   Input Parameters:
270: + da  - The `DMDA`
271: . xo  - The offset in the x direction
272: . yo  - The offset in the y direction
273: - zo  - The offset in the z direction

275:   Level: intermediate

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

281: .seealso: `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
282: @*/
283: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
284: {
285:   DM_DA *dd = (DM_DA *)da->data;

294:   dd->xo = xo;
295:   dd->yo = yo;
296:   dd->zo = zo;
297:   dd->Mo = Mo;
298:   dd->No = No;
299:   dd->Po = Po;

301:   if (da->coordinates[0].dm) DMDASetOffset(da->coordinates[0].dm, xo, yo, zo, Mo, No, Po);
302:   return 0;
303: }

305: /*@
306:   DMDAGetOffset - Gets the index offset of the `DMDA`.

308:   Not collective

310:   Input Parameter:
311: . da  - The `DMDA`

313:   Output Parameters:
314: + xo  - The offset in the x direction
315: . yo  - The offset in the y direction
316: . zo  - The offset in the z direction
317: . Mo  - The global size in the x direction
318: . No  - The global size in the y direction
319: - Po  - The global size in the z direction

321:   Level: intermediate

323: .seealso: `DM`, `DMDA`, `DMDASetOffset()`, `DMDAVecGetArray()`
324: @*/
325: PetscErrorCode DMDAGetOffset(DM da, PetscInt *xo, PetscInt *yo, PetscInt *zo, PetscInt *Mo, PetscInt *No, PetscInt *Po)
326: {
327:   DM_DA *dd = (DM_DA *)da->data;

330:   if (xo) *xo = dd->xo;
331:   if (yo) *yo = dd->yo;
332:   if (zo) *zo = dd->zo;
333:   if (Mo) *Mo = dd->Mo;
334:   if (No) *No = dd->No;
335:   if (Po) *Po = dd->Po;
336:   return 0;
337: }

339: /*@
340:   DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain `DM`.

342:   Not collective

344:   Input Parameter:
345: . da  - The `DMDA`

347:   Output Parameters:
348: + xs  - The start of the region in x
349: . ys  - The start of the region in y
350: . zs  - The start of the region in z
351: . xs  - The size of the region in x
352: . ys  - The size of the region in y
353: - zs  - The size of the region in z

355:   Level: intermediate

357: .seealso: `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
358: @*/
359: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
360: {
361:   DM_DA *dd = (DM_DA *)da->data;

364:   if (xs) *xs = dd->nonxs;
365:   if (ys) *ys = dd->nonys;
366:   if (zs) *zs = dd->nonzs;
367:   if (xm) *xm = dd->nonxm;
368:   if (ym) *ym = dd->nonym;
369:   if (zm) *zm = dd->nonzm;
370:   return 0;
371: }

373: /*@
374:   DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain `DM`.

376:   Collective on da

378:   Input Parameters:
379: + da  - The `DMDA`
380: . xs  - The start of the region in x
381: . ys  - The start of the region in y
382: . zs  - The start of the region in z
383: . xs  - The size of the region in x
384: . ys  - The size of the region in y
385: - zs  - The size of the region in z

387:   Level: intermediate

389: .seealso: `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
390: @*/
391: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
392: {
393:   DM_DA *dd = (DM_DA *)da->data;

402:   dd->nonxs = xs;
403:   dd->nonys = ys;
404:   dd->nonzs = zs;
405:   dd->nonxm = xm;
406:   dd->nonym = ym;
407:   dd->nonzm = zm;

409:   return 0;
410: }

412: /*@
413:   DMDASetStencilType - Sets the type of the communication stencil

415:   Logically Collective on da

417:   Input Parameters:
418: + da    - The `DMDA`
419: - stype - The stencil type, use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.

421:   Level: intermediate

423: .seealso: `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`
424: @*/
425: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
426: {
427:   DM_DA *dd = (DM_DA *)da->data;

432:   dd->stencil_type = stype;
433:   return 0;
434: }

436: /*@
437:   DMDAGetStencilType - Gets the type of the communication stencil

439:   Not collective

441:   Input Parameter:
442: . da    - The `DMDA`

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

447:   Level: intermediate

449: .seealso: `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`
450: @*/
451: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
452: {
453:   DM_DA *dd = (DM_DA *)da->data;

457:   *stype = dd->stencil_type;
458:   return 0;
459: }

461: /*@
462:   DMDASetStencilWidth - Sets the width of the communication stencil

464:   Logically Collective on da

466:   Input Parameters:
467: + da    - The `DMDA`
468: - width - The stencil width

470:   Level: intermediate

472: .seealso: `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`
473: @*/
474: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
475: {
476:   DM_DA *dd = (DM_DA *)da->data;

481:   dd->s = width;
482:   return 0;
483: }

485: /*@
486:   DMDAGetStencilWidth - Gets the width of the communication stencil

488:   Not collective

490:   Input Parameter:
491: . da    - The `DMDA`

493:   Output Parameter:
494: . width - The stencil width

496:   Level: intermediate

498: .seealso: `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`
499: @*/
500: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
501: {
502:   DM_DA *dd = (DM_DA *)da->data;

506:   *width = dd->s;
507:   return 0;
508: }

510: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da, PetscInt M, PetscInt m, const PetscInt lx[])
511: {
512:   PetscInt i, sum;

515:   for (i = sum = 0; i < m; i++) sum += lx[i];
517:   return 0;
518: }

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

523:   Logically Collective on da

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

531:   Level: intermediate

533:   Note: these numbers are NOT multiplied by the number of dof per node.

535: .seealso: `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`
536: @*/
537: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
538: {
539:   DM_DA *dd = (DM_DA *)da->data;

543:   if (lx) {
545:     DMDACheckOwnershipRanges_Private(da, dd->M, dd->m, lx);
546:     if (!dd->lx) PetscMalloc1(dd->m, &dd->lx);
547:     PetscArraycpy(dd->lx, lx, dd->m);
548:   }
549:   if (ly) {
551:     DMDACheckOwnershipRanges_Private(da, dd->N, dd->n, ly);
552:     if (!dd->ly) PetscMalloc1(dd->n, &dd->ly);
553:     PetscArraycpy(dd->ly, ly, dd->n);
554:   }
555:   if (lz) {
557:     DMDACheckOwnershipRanges_Private(da, dd->P, dd->p, lz);
558:     if (!dd->lz) PetscMalloc1(dd->p, &dd->lz);
559:     PetscArraycpy(dd->lz, lz, dd->p);
560:   }
561:   return 0;
562: }

564: /*@
565:        DMDASetInterpolationType - Sets the type of interpolation that will be
566:           returned by `DMCreateInterpolation()`

568:    Logically Collective on da

570:    Input Parameters:
571: +  da - initial distributed array
572: -  ctype - `DMDA_Q1` and `DMDA_Q0` are currently the only supported forms

574:    Level: intermediate

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

579: .seealso: `DM`, `DMDA`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMDestroy()`, `DMDAInterpolationType`
580: @*/
581: PetscErrorCode DMDASetInterpolationType(DM da, DMDAInterpolationType ctype)
582: {
583:   DM_DA *dd = (DM_DA *)da->data;

587:   dd->interptype = ctype;
588:   return 0;
589: }

591: /*@
592:        DMDAGetInterpolationType - Gets the type of interpolation that will be
593:           used by `DMCreateInterpolation()`

595:    Not Collective

597:    Input Parameter:
598: .  da - distributed array

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

603:    Level: intermediate

605: .seealso: `DM`, `DMDA`, `DMDAInterpolationType`, `DMDASetInterpolationType()`, `DMCreateInterpolation()`
606: @*/
607: PetscErrorCode DMDAGetInterpolationType(DM da, DMDAInterpolationType *ctype)
608: {
609:   DM_DA *dd = (DM_DA *)da->data;

613:   *ctype = dd->interptype;
614:   return 0;
615: }

617: /*@C
618:       DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
619:         processes neighbors.

621:     Not Collective

623:    Input Parameter:
624: .     da - the `DMDA` object

626:    Output Parameters:
627: .     ranks - the neighbors ranks, stored with the x index increasing most rapidly.
628:               this process itself is in the list

630:    Level: intermediate

632:    Notes:
633:    In 2d the array is of length 9, in 3d of length 27

635:    Not supported in 1d

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

639:    Fortran Note:
640:     In fortran you must pass in an array of the appropriate length.

642: .seealso: `DMDA`, `DM`
643: @*/
644: PetscErrorCode DMDAGetNeighbors(DM da, const PetscMPIInt *ranks[])
645: {
646:   DM_DA *dd = (DM_DA *)da->data;

649:   *ranks = dd->neighbors;
650:   return 0;
651: }

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

656:     Not Collective

658:    Input Parameter:
659: .     da - the `DMDA` object

661:    Output Parameters:
662: +     lx - ownership along x direction (optional)
663: .     ly - ownership along y direction (optional)
664: -     lz - ownership along z direction (optional)

666:    Level: intermediate

668:     Note:
669:     These correspond to the optional final arguments passed to `DMDACreate()`, `DMDACreate2d()`, `DMDACreate3d()`

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

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

676:     Fortran Note:
677:     In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
678:     eighth arguments from `DMDAGetInfo()`

680: .seealso: `DM`, `DMDA`, `DMDAGetCorners()`, `DMDAGetGhostCorners()`, `DMDACreate()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `VecGetOwnershipRanges()`
681: @*/
682: PetscErrorCode DMDAGetOwnershipRanges(DM da, const PetscInt *lx[], const PetscInt *ly[], const PetscInt *lz[])
683: {
684:   DM_DA *dd = (DM_DA *)da->data;

687:   if (lx) *lx = dd->lx;
688:   if (ly) *ly = dd->ly;
689:   if (lz) *lz = dd->lz;
690:   return 0;
691: }

693: /*@
694:      DMDASetRefinementFactor - Set the ratios that the `DMDA` grid is refined

696:     Logically Collective on da

698:   Input Parameters:
699: +    da - the `DMDA` object
700: .    refine_x - ratio of fine grid to coarse in x direction (2 by default)
701: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
702: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

704:   Options Database Keys:
705: +  -da_refine_x refine_x - refinement ratio in x direction
706: .  -da_refine_y rafine_y - refinement ratio in y direction
707: .  -da_refine_z refine_z - refinement ratio in z direction
708: -  -da_refine <n> - refine the DMDA object n times when it is created.

710:   Level: intermediate

712:     Note:
713:     Pass `PETSC_IGNORE` to leave a value unchanged

715: .seealso: `DM`, `DMDA`, `DMRefine()`, `DMDAGetRefinementFactor()`
716: @*/
717: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y, PetscInt refine_z)
718: {
719:   DM_DA *dd = (DM_DA *)da->data;


726:   if (refine_x > 0) dd->refine_x = refine_x;
727:   if (refine_y > 0) dd->refine_y = refine_y;
728:   if (refine_z > 0) dd->refine_z = refine_z;
729:   return 0;
730: }

732: /*@C
733:      DMDAGetRefinementFactor - Gets the ratios that the `DMDA` grid is refined

735:     Not Collective

737:   Input Parameter:
738: .    da - the `DMDA` object

740:   Output Parameters:
741: +    refine_x - ratio of fine grid to coarse in x direction (2 by default)
742: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
743: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

745:   Level: intermediate

747:     Note:
748:     Pass NULL for values you do not need

750: .seealso: `DM`, `DMDA`, `DMRefine()`, `DMDASetRefinementFactor()`
751: @*/
752: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y, PetscInt *refine_z)
753: {
754:   DM_DA *dd = (DM_DA *)da->data;

757:   if (refine_x) *refine_x = dd->refine_x;
758:   if (refine_y) *refine_y = dd->refine_y;
759:   if (refine_z) *refine_z = dd->refine_z;
760:   return 0;
761: }

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

766:     Logically Collective on da

768:   Input Parameters:
769: +    da - the `DMDA` object
770: -    f - the function that allocates the matrix for that specific DMDA

772:   Level: developer

774:    Note:
775:     See `DMDASetBlockFills()` that provides a simple way to provide the nonzero structure for
776:        the diagonal and off-diagonal blocks of the matrix

778:    Fortran Note:
779:    Not supported from Fortran

781: .seealso: `DM`, `DMDA`, `DMCreateMatrix()`, `DMDASetBlockFills()`
782: @*/
783: PetscErrorCode DMDASetGetMatrix(DM da, PetscErrorCode (*f)(DM, Mat *))
784: {
786:   da->ops->creatematrix = f;
787:   return 0;
788: }

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

793:    Not Collective

795:    Input Parameters:
796: +  da - the `DMDA` object
797: .  m - number of MatStencils
798: -  idxm - grid points (and component number when dof > 1)

800:    Output Parameter:
801: .  gidxm - global row indices

803:    Level: intermediate

805: .seealso: `DM`, `DMDA`, `MatStencil`
806: @*/
807: PetscErrorCode DMDAMapMatStencilToGlobal(DM da, PetscInt m, const MatStencil idxm[], PetscInt gidxm[])
808: {
809:   const DM_DA           *dd  = (const DM_DA *)da->data;
810:   const PetscInt        *dxm = (const PetscInt *)idxm;
811:   PetscInt               i, j, sdim, tmp, dim;
812:   PetscInt               dims[4], starts[4], dims2[3], starts2[3], dof = dd->w;
813:   ISLocalToGlobalMapping ltog;

815:   if (m <= 0) return 0;

817:   /* Code adapted from DMDAGetGhostCorners() */
818:   starts2[0] = dd->Xs / dof + dd->xo;
819:   starts2[1] = dd->Ys + dd->yo;
820:   starts2[2] = dd->Zs + dd->zo;
821:   dims2[0]   = (dd->Xe - dd->Xs) / dof;
822:   dims2[1]   = (dd->Ye - dd->Ys);
823:   dims2[2]   = (dd->Ze - dd->Zs);

825:   /* As if we do MatSetStencil() to get dims[]/starts[] of mat->stencil */
826:   dim  = da->dim;                   /* DA dim: 1 to 3 */
827:   sdim = dim + (dof > 1 ? 1 : 0);   /* Dimensions in MatStencil's (k,j,i,c) view */
828:   for (i = 0; i < dim; i++) {       /* Reverse the order and also skip the unused dimensions */
829:     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 */
830:     starts[i] = starts2[dim - i - 1];
831:   }
832:   starts[dim] = 0; /* Append the extra dim for dof (won't be used below if dof=1) */
833:   dims[dim]   = dof;

835:   /* Map stencils to local indices (code adapted from MatSetValuesStencil()) */
836:   for (i = 0; i < m; i++) {
837:     dxm += 3 - dim; /* Input is {k,j,i,c}; move the pointer to the first used index, e.g., j in 2D */
838:     tmp = 0;
839:     for (j = 0; j < sdim; j++) {                                                      /* Iter over, ex. j,i or j,i,c in 2D */
840:       if (tmp < 0 || dxm[j] < starts[j] || dxm[j] >= (starts[j] + dims[j])) tmp = -1; /* Beyond the ghost region, therefore ignored with negative indices */
841:       else tmp = tmp * dims[j] + (dxm[j] - starts[j]);
842:     }
843:     gidxm[i] = tmp;
844:     /* Move to the next MatStencil point */
845:     if (dof > 1) dxm += sdim; /* c is already counted in sdim */
846:     else dxm += sdim + 1;     /* skip the unused c */
847:   }

849:   /* Map local indices to global indices */
850:   DMGetLocalToGlobalMapping(da, &ltog);
851:   ISLocalToGlobalMappingApply(ltog, m, gidxm, gidxm);
852:   return 0;
853: }

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

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

866:   if (ratio == 1) {
867:     PetscArraycpy(lf, lc, m);
868:     return 0;
869:   }
870:   for (i = 0; i < m; i++) totalc += lc[i];
871:   remaining = (!periodic) + ratio * (totalc - (!periodic));
872:   for (i = 0; i < m; i++) {
873:     PetscInt want = remaining / (m - i) + !!(remaining % (m - i));
874:     if (i == m - 1) lf[i] = want;
875:     else {
876:       const PetscInt nextc = startc + lc[i];
877:       /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
878:        * coarse stencil width of the first coarse node in the next subdomain. */
879:       while ((startf + want) / ratio < nextc - stencil_width) want++;
880:       /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
881:        * coarse stencil width of the last coarse node in the current subdomain. */
882:       while ((startf + want - 1 + ratio - 1) / ratio > nextc - 1 + stencil_width) want--;
883:       /* Make sure all constraints are satisfied */
884:       if (want < 0 || want > remaining || ((startf + want) / ratio < nextc - stencil_width) || ((startf + want - 1 + ratio - 1) / ratio > nextc - 1 + stencil_width))
885:         SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Could not find a compatible refined ownership range");
886:     }
887:     lf[i] = want;
888:     startc += lc[i];
889:     startf += lf[i];
890:     remaining -= lf[i];
891:   }
892:   return 0;
893: }

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

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

906:   if (ratio == 1) {
907:     PetscArraycpy(lc, lf, m);
908:     return 0;
909:   }
910:   for (i = 0, totalf = 0; i < m; i++) totalf += lf[i];
911:   remaining = (!periodic) + (totalf - (!periodic)) / ratio;
912:   for (i = 0, startc = 0, startf = 0; i < m; i++) {
913:     PetscInt want = remaining / (m - i) + !!(remaining % (m - i));
914:     if (i == m - 1) lc[i] = want;
915:     else {
916:       const PetscInt nextf = startf + lf[i];
917:       /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
918:        * node is within one stencil width. */
919:       while (nextf / ratio < startc + want - stencil_width) want--;
920:       /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
921:        * fine node is within one stencil width. */
922:       while ((nextf - 1 + ratio - 1) / ratio > startc + want - 1 + stencil_width) want++;
923:       if (want < 0 || want > remaining || (nextf / ratio < startc + want - stencil_width) || ((nextf - 1 + ratio - 1) / ratio > startc + want - 1 + stencil_width))
924:         SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Could not find a compatible coarsened ownership range");
925:     }
926:     lc[i] = want;
927:     startc += lc[i];
928:     startf += lf[i];
929:     remaining -= lc[i];
930:   }
931:   return 0;
932: }

934: PetscErrorCode DMRefine_DA(DM da, MPI_Comm comm, DM *daref)
935: {
936:   PetscInt M, N, P, i, dim;
937:   Vec      coordsc, coordsf;
938:   DM       da2;
939:   DM_DA   *dd = (DM_DA *)da->data, *dd2;


944:   DMGetDimension(da, &dim);
945:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
946:     M = dd->refine_x * dd->M;
947:   } else {
948:     M = 1 + dd->refine_x * (dd->M - 1);
949:   }
950:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
951:     if (dim > 1) {
952:       N = dd->refine_y * dd->N;
953:     } else {
954:       N = 1;
955:     }
956:   } else {
957:     N = 1 + dd->refine_y * (dd->N - 1);
958:   }
959:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
960:     if (dim > 2) {
961:       P = dd->refine_z * dd->P;
962:     } else {
963:       P = 1;
964:     }
965:   } else {
966:     P = 1 + dd->refine_z * (dd->P - 1);
967:   }
968:   DMDACreate(PetscObjectComm((PetscObject)da), &da2);
969:   DMSetOptionsPrefix(da2, ((PetscObject)da)->prefix);
970:   DMSetDimension(da2, dim);
971:   DMDASetSizes(da2, M, N, P);
972:   DMDASetNumProcs(da2, dd->m, dd->n, dd->p);
973:   DMDASetBoundaryType(da2, dd->bx, dd->by, dd->bz);
974:   DMDASetDof(da2, dd->w);
975:   DMDASetStencilType(da2, dd->stencil_type);
976:   DMDASetStencilWidth(da2, dd->s);
977:   if (dim == 3) {
978:     PetscInt *lx, *ly, *lz;
979:     PetscMalloc3(dd->m, &lx, dd->n, &ly, dd->p, &lz);
980:     DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx);
981:     DMDARefineOwnershipRanges(da, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_y, dd->n, dd->ly, ly);
982:     DMDARefineOwnershipRanges(da, (PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_z, dd->p, dd->lz, lz);
983:     DMDASetOwnershipRanges(da2, lx, ly, lz);
984:     PetscFree3(lx, ly, lz);
985:   } else if (dim == 2) {
986:     PetscInt *lx, *ly;
987:     PetscMalloc2(dd->m, &lx, dd->n, &ly);
988:     DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx);
989:     DMDARefineOwnershipRanges(da, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_y, dd->n, dd->ly, ly);
990:     DMDASetOwnershipRanges(da2, lx, ly, NULL);
991:     PetscFree2(lx, ly);
992:   } else if (dim == 1) {
993:     PetscInt *lx;
994:     PetscMalloc1(dd->m, &lx);
995:     DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx);
996:     DMDASetOwnershipRanges(da2, lx, NULL, NULL);
997:     PetscFree(lx);
998:   }
999:   dd2 = (DM_DA *)da2->data;

1001:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
1002:   da2->ops->creatematrix = da->ops->creatematrix;
1003:   /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
1004:   da2->ops->getcoloring = da->ops->getcoloring;
1005:   dd2->interptype       = dd->interptype;

1007:   /* copy fill information if given */
1008:   if (dd->dfill) {
1009:     PetscMalloc1(dd->dfill[dd->w] + dd->w + 1, &dd2->dfill);
1010:     PetscArraycpy(dd2->dfill, dd->dfill, dd->dfill[dd->w] + dd->w + 1);
1011:   }
1012:   if (dd->ofill) {
1013:     PetscMalloc1(dd->ofill[dd->w] + dd->w + 1, &dd2->ofill);
1014:     PetscArraycpy(dd2->ofill, dd->ofill, dd->ofill[dd->w] + dd->w + 1);
1015:   }
1016:   /* copy the refine information */
1017:   dd2->coarsen_x = dd2->refine_x = dd->refine_x;
1018:   dd2->coarsen_y = dd2->refine_y = dd->refine_y;
1019:   dd2->coarsen_z = dd2->refine_z = dd->refine_z;

1021:   if (dd->refine_z_hier) {
1022:     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];
1023:     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];
1024:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1025:     PetscMalloc1(dd2->refine_z_hier_n, &dd2->refine_z_hier);
1026:     PetscArraycpy(dd2->refine_z_hier, dd->refine_z_hier, dd2->refine_z_hier_n);
1027:   }
1028:   if (dd->refine_y_hier) {
1029:     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];
1030:     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];
1031:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1032:     PetscMalloc1(dd2->refine_y_hier_n, &dd2->refine_y_hier);
1033:     PetscArraycpy(dd2->refine_y_hier, dd->refine_y_hier, dd2->refine_y_hier_n);
1034:   }
1035:   if (dd->refine_x_hier) {
1036:     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];
1037:     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];
1038:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1039:     PetscMalloc1(dd2->refine_x_hier_n, &dd2->refine_x_hier);
1040:     PetscArraycpy(dd2->refine_x_hier, dd->refine_x_hier, dd2->refine_x_hier_n);
1041:   }

1043:   /* copy vector type information */
1044:   DMSetVecType(da2, da->vectype);

1046:   dd2->lf = dd->lf;
1047:   dd2->lj = dd->lj;

1049:   da2->leveldown = da->leveldown;
1050:   da2->levelup   = da->levelup + 1;

1052:   DMSetUp(da2);

1054:   /* interpolate coordinates if they are set on the coarse grid */
1055:   DMGetCoordinates(da, &coordsc);
1056:   if (coordsc) {
1057:     DM  cdaf, cdac;
1058:     Mat II;

1060:     DMGetCoordinateDM(da, &cdac);
1061:     DMGetCoordinateDM(da2, &cdaf);
1062:     /* force creation of the coordinate vector */
1063:     DMDASetUniformCoordinates(da2, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
1064:     DMGetCoordinates(da2, &coordsf);
1065:     DMCreateInterpolation(cdac, cdaf, &II, NULL);
1066:     MatInterpolate(II, coordsc, coordsf);
1067:     MatDestroy(&II);
1068:   }

1070:   for (i = 0; i < da->bs; i++) {
1071:     const char *fieldname;
1072:     DMDAGetFieldName(da, i, &fieldname);
1073:     DMDASetFieldName(da2, i, fieldname);
1074:   }

1076:   *daref = da2;
1077:   return 0;
1078: }

1080: PetscErrorCode DMCoarsen_DA(DM dmf, MPI_Comm comm, DM *dmc)
1081: {
1082:   PetscInt M, N, P, i, dim;
1083:   Vec      coordsc, coordsf;
1084:   DM       dmc2;
1085:   DM_DA   *dd = (DM_DA *)dmf->data, *dd2;


1090:   DMGetDimension(dmf, &dim);
1091:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1092:     M = dd->M / dd->coarsen_x;
1093:   } else {
1094:     M = 1 + (dd->M - 1) / dd->coarsen_x;
1095:   }
1096:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1097:     if (dim > 1) {
1098:       N = dd->N / dd->coarsen_y;
1099:     } else {
1100:       N = 1;
1101:     }
1102:   } else {
1103:     N = 1 + (dd->N - 1) / dd->coarsen_y;
1104:   }
1105:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1106:     if (dim > 2) {
1107:       P = dd->P / dd->coarsen_z;
1108:     } else {
1109:       P = 1;
1110:     }
1111:   } else {
1112:     P = 1 + (dd->P - 1) / dd->coarsen_z;
1113:   }
1114:   DMDACreate(PetscObjectComm((PetscObject)dmf), &dmc2);
1115:   DMSetOptionsPrefix(dmc2, ((PetscObject)dmf)->prefix);
1116:   DMSetDimension(dmc2, dim);
1117:   DMDASetSizes(dmc2, M, N, P);
1118:   DMDASetNumProcs(dmc2, dd->m, dd->n, dd->p);
1119:   DMDASetBoundaryType(dmc2, dd->bx, dd->by, dd->bz);
1120:   DMDASetDof(dmc2, dd->w);
1121:   DMDASetStencilType(dmc2, dd->stencil_type);
1122:   DMDASetStencilWidth(dmc2, dd->s);
1123:   if (dim == 3) {
1124:     PetscInt *lx, *ly, *lz;
1125:     PetscMalloc3(dd->m, &lx, dd->n, &ly, dd->p, &lz);
1126:     DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx);
1127:     DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_y, dd->n, dd->ly, ly);
1128:     DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_z, dd->p, dd->lz, lz);
1129:     DMDASetOwnershipRanges(dmc2, lx, ly, lz);
1130:     PetscFree3(lx, ly, lz);
1131:   } else if (dim == 2) {
1132:     PetscInt *lx, *ly;
1133:     PetscMalloc2(dd->m, &lx, dd->n, &ly);
1134:     DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx);
1135:     DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_y, dd->n, dd->ly, ly);
1136:     DMDASetOwnershipRanges(dmc2, lx, ly, NULL);
1137:     PetscFree2(lx, ly);
1138:   } else if (dim == 1) {
1139:     PetscInt *lx;
1140:     PetscMalloc1(dd->m, &lx);
1141:     DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx);
1142:     DMDASetOwnershipRanges(dmc2, lx, NULL, NULL);
1143:     PetscFree(lx);
1144:   }
1145:   dd2 = (DM_DA *)dmc2->data;

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

1153:   /* copy fill information if given */
1154:   if (dd->dfill) {
1155:     PetscMalloc1(dd->dfill[dd->w] + dd->w + 1, &dd2->dfill);
1156:     PetscArraycpy(dd2->dfill, dd->dfill, dd->dfill[dd->w] + dd->w + 1);
1157:   }
1158:   if (dd->ofill) {
1159:     PetscMalloc1(dd->ofill[dd->w] + dd->w + 1, &dd2->ofill);
1160:     PetscArraycpy(dd2->ofill, dd->ofill, dd->ofill[dd->w] + dd->w + 1);
1161:   }
1162:   /* copy the refine information */
1163:   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1164:   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1165:   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;

1167:   if (dd->refine_z_hier) {
1168:     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];
1169:     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];
1170:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1171:     PetscMalloc1(dd2->refine_z_hier_n, &dd2->refine_z_hier);
1172:     PetscArraycpy(dd2->refine_z_hier, dd->refine_z_hier, dd2->refine_z_hier_n);
1173:   }
1174:   if (dd->refine_y_hier) {
1175:     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];
1176:     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];
1177:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1178:     PetscMalloc1(dd2->refine_y_hier_n, &dd2->refine_y_hier);
1179:     PetscArraycpy(dd2->refine_y_hier, dd->refine_y_hier, dd2->refine_y_hier_n);
1180:   }
1181:   if (dd->refine_x_hier) {
1182:     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];
1183:     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];
1184:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1185:     PetscMalloc1(dd2->refine_x_hier_n, &dd2->refine_x_hier);
1186:     PetscArraycpy(dd2->refine_x_hier, dd->refine_x_hier, dd2->refine_x_hier_n);
1187:   }

1189:   /* copy vector type information */
1190:   DMSetVecType(dmc2, dmf->vectype);

1192:   dd2->lf = dd->lf;
1193:   dd2->lj = dd->lj;

1195:   dmc2->leveldown = dmf->leveldown + 1;
1196:   dmc2->levelup   = dmf->levelup;

1198:   DMSetUp(dmc2);

1200:   /* inject coordinates if they are set on the fine grid */
1201:   DMGetCoordinates(dmf, &coordsf);
1202:   if (coordsf) {
1203:     DM         cdaf, cdac;
1204:     Mat        inject;
1205:     VecScatter vscat;

1207:     DMGetCoordinateDM(dmf, &cdaf);
1208:     DMGetCoordinateDM(dmc2, &cdac);
1209:     /* force creation of the coordinate vector */
1210:     DMDASetUniformCoordinates(dmc2, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
1211:     DMGetCoordinates(dmc2, &coordsc);

1213:     DMCreateInjection(cdac, cdaf, &inject);
1214:     MatScatterGetVecScatter(inject, &vscat);
1215:     VecScatterBegin(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD);
1216:     VecScatterEnd(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD);
1217:     MatDestroy(&inject);
1218:   }

1220:   for (i = 0; i < dmf->bs; i++) {
1221:     const char *fieldname;
1222:     DMDAGetFieldName(dmf, i, &fieldname);
1223:     DMDASetFieldName(dmc2, i, fieldname);
1224:   }

1226:   *dmc = dmc2;
1227:   return 0;
1228: }

1230: PetscErrorCode DMRefineHierarchy_DA(DM da, PetscInt nlevels, DM daf[])
1231: {
1232:   PetscInt i, n, *refx, *refy, *refz;

1236:   if (nlevels == 0) return 0;

1239:   /* Get refinement factors, defaults taken from the coarse DMDA */
1240:   PetscMalloc3(nlevels, &refx, nlevels, &refy, nlevels, &refz);
1241:   for (i = 0; i < nlevels; i++) DMDAGetRefinementFactor(da, &refx[i], &refy[i], &refz[i]);
1242:   n = nlevels;
1243:   PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_x", refx, &n, NULL);
1244:   n = nlevels;
1245:   PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_y", refy, &n, NULL);
1246:   n = nlevels;
1247:   PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_z", refz, &n, NULL);

1249:   DMDASetRefinementFactor(da, refx[0], refy[0], refz[0]);
1250:   DMRefine(da, PetscObjectComm((PetscObject)da), &daf[0]);
1251:   for (i = 1; i < nlevels; i++) {
1252:     DMDASetRefinementFactor(daf[i - 1], refx[i], refy[i], refz[i]);
1253:     DMRefine(daf[i - 1], PetscObjectComm((PetscObject)da), &daf[i]);
1254:   }
1255:   PetscFree3(refx, refy, refz);
1256:   return 0;
1257: }

1259: PetscErrorCode DMCoarsenHierarchy_DA(DM da, PetscInt nlevels, DM dac[])
1260: {
1261:   PetscInt i;

1265:   if (nlevels == 0) return 0;
1267:   DMCoarsen(da, PetscObjectComm((PetscObject)da), &dac[0]);
1268:   for (i = 1; i < nlevels; i++) DMCoarsen(dac[i - 1], PetscObjectComm((PetscObject)da), &dac[i]);
1269:   return 0;
1270: }

1272: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm, PetscInt n, PetscReal *nodes)
1273: {
1274:   PetscInt     i, j, xs, xn, q;
1275:   PetscScalar *xx;
1276:   PetscReal    h;
1277:   Vec          x;
1278:   DM_DA       *da = (DM_DA *)dm->data;

1280:   if (da->bx != DM_BOUNDARY_PERIODIC) {
1281:     DMDAGetInfo(dm, NULL, &q, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
1282:     q = (q - 1) / (n - 1); /* number of spectral elements */
1283:     h = 2.0 / q;
1284:     DMDAGetCorners(dm, &xs, NULL, NULL, &xn, NULL, NULL);
1285:     xs = xs / (n - 1);
1286:     xn = xn / (n - 1);
1287:     DMDASetUniformCoordinates(dm, -1., 1., 0., 0., 0., 0.);
1288:     DMGetCoordinates(dm, &x);
1289:     DMDAVecGetArray(dm, x, &xx);

1291:     /* loop over local spectral elements */
1292:     for (j = xs; j < xs + xn; j++) {
1293:       /*
1294:        Except for the first process, each process starts on the second GLL point of the first element on that process
1295:        */
1296:       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.;
1297:     }
1298:     DMDAVecRestoreArray(dm, x, &xx);
1299:   } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not yet implemented for periodic");
1300:   return 0;
1301: }

1303: /*@

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

1307:    Collective on da

1309:    Input Parameters:
1310: +   da - the `DMDA` object
1311: -   n - the number of GLL nodes
1312: -   nodes - the GLL nodes

1314:    Level: advanced

1316:    Note:
1317:    The parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1318:    on each process much be divisible by the number of GLL elements needed per process. This depends on whether the `DM` is
1319:    periodic or not.

1321: .seealso: `DM`, `DMDA`, `DMDACreate()`, `PetscDTGaussLobattoLegendreQuadrature()`, `DMGetCoordinates()`
1322: @*/
1323: PetscErrorCode DMDASetGLLCoordinates(DM da, PetscInt n, PetscReal *nodes)
1324: {
1325:   if (da->dim == 1) {
1326:     DMDASetGLLCoordinates_1d(da, n, nodes);
1327:   } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not yet implemented for 2 or 3d");
1328:   return 0;
1329: }

1331: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM da1, DM dm2, PetscBool *compatible, PetscBool *set)
1332: {
1333:   DM_DA    *dd1 = (DM_DA *)da1->data, *dd2;
1334:   DM        da2;
1335:   DMType    dmtype2;
1336:   PetscBool isda, compatibleLocal;
1337:   PetscInt  i;

1340:   DMGetType(dm2, &dmtype2);
1341:   PetscStrcmp(dmtype2, DMDA, &isda);
1342:   if (isda) {
1343:     da2 = dm2;
1344:     dd2 = (DM_DA *)da2->data;
1346:     compatibleLocal = (PetscBool)(da1->dim == da2->dim);
1347:     if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
1348:     /*                                                                           Global size              ranks               Boundary type */
1349:     if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
1350:     if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
1351:     if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
1352:     if (compatibleLocal) {
1353:       for (i = 0; i < dd1->m; ++i) { compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i])); /* Local size     */ }
1354:     }
1355:     if (compatibleLocal && da1->dim > 1) {
1356:       for (i = 0; i < dd1->n; ++i) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
1357:     }
1358:     if (compatibleLocal && da1->dim > 2) {
1359:       for (i = 0; i < dd1->p; ++i) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
1360:     }
1361:     *compatible = compatibleLocal;
1362:     *set        = PETSC_TRUE;
1363:   } else {
1364:     /* Decline to determine compatibility with other DM types */
1365:     *set = PETSC_FALSE;
1366:   }
1367:   return 0;
1368: }