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 Notes:
 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: `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: `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: `DMDACreate()`, `DMDestroy()`, `DMDA`, `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: `DMDAGetDof()`, `DMDACreate()`, `DMDestroy()`, `DMDA`
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: `DMDASetDof()`, `DMDACreate()`, `DMDestroy()`, `DMDA`
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: `DMCreateDomainDecomposition()`, `DMDASetOverlap()`, `DMDA`
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: `DMCreateDomainDecomposition()`, `DMDAGetOverlap()`, `DMDA`
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: `DMCreateDomainDecomposition()`, `DMDASetNumLocalSubDomains()`, `DMDA`
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: `DMCreateDomainDecomposition()`, `DMDAGetNumLocalSubDomains()`, `DMDA`
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:   Notes:
278:     This is used primarily to overlap a computation on a local DA with that on a global DA without
279:   changing boundary conditions or subdomain features that depend upon the global offsets.

281: .seealso: `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 DA.

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: `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: `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: `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: `DMDACreate()`, `DMDestroy()`, `DMDA`
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: `DMDACreate()`, `DMDestroy()`, `DMDA`
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: `DMDACreate()`, `DMDestroy()`, `DMDA`
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: `DMDACreate()`, `DMDestroy()`, `DMDA`
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: `DMDACreate()`, `DMDestroy()`, `DMDA`
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:    Notes:
577:     you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()

579: .seealso: `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMDestroy()`, `DMDA`, `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: `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:    Notes:
631:     In 2d the array is of length 9, in 3d of length 27
632:           Not supported in 1d
633:           Do not free the array, it is freed when the DMDA is destroyed.

635:    Fortran Notes:
636:     In fortran you must pass in an array of the appropriate length.

638:    Level: intermediate

640: @*/
641: PetscErrorCode DMDAGetNeighbors(DM da, const PetscMPIInt *ranks[])
642: {
643:   DM_DA *dd = (DM_DA *)da->data;

646:   *ranks = dd->neighbors;
647:   return 0;
648: }

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

653:     Not Collective

655:    Input Parameter:
656: .     da - the DMDA object

658:    Output Parameters:
659: +     lx - ownership along x direction (optional)
660: .     ly - ownership along y direction (optional)
661: -     lz - ownership along z direction (optional)

663:    Level: intermediate

665:     Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()

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

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

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

675: .seealso: `DMDAGetCorners()`, `DMDAGetGhostCorners()`, `DMDACreate()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `VecGetOwnershipRanges()`
676: @*/
677: PetscErrorCode DMDAGetOwnershipRanges(DM da, const PetscInt *lx[], const PetscInt *ly[], const PetscInt *lz[])
678: {
679:   DM_DA *dd = (DM_DA *)da->data;

682:   if (lx) *lx = dd->lx;
683:   if (ly) *ly = dd->ly;
684:   if (lz) *lz = dd->lz;
685:   return 0;
686: }

688: /*@
689:      DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined

691:     Logically Collective on da

693:   Input Parameters:
694: +    da - the DMDA object
695: .    refine_x - ratio of fine grid to coarse in x direction (2 by default)
696: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
697: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

699:   Options Database:
700: +  -da_refine_x refine_x - refinement ratio in x direction
701: .  -da_refine_y rafine_y - refinement ratio in y direction
702: .  -da_refine_z refine_z - refinement ratio in z direction
703: -  -da_refine <n> - refine the DMDA object n times when it is created.

705:   Level: intermediate

707:     Notes:
708:     Pass PETSC_IGNORE to leave a value unchanged

710: .seealso: `DMRefine()`, `DMDAGetRefinementFactor()`
711: @*/
712: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y, PetscInt refine_z)
713: {
714:   DM_DA *dd = (DM_DA *)da->data;


721:   if (refine_x > 0) dd->refine_x = refine_x;
722:   if (refine_y > 0) dd->refine_y = refine_y;
723:   if (refine_z > 0) dd->refine_z = refine_z;
724:   return 0;
725: }

727: /*@C
728:      DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined

730:     Not Collective

732:   Input Parameter:
733: .    da - the DMDA object

735:   Output Parameters:
736: +    refine_x - ratio of fine grid to coarse in x direction (2 by default)
737: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
738: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

740:   Level: intermediate

742:     Notes:
743:     Pass NULL for values you do not need

745: .seealso: `DMRefine()`, `DMDASetRefinementFactor()`
746: @*/
747: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y, PetscInt *refine_z)
748: {
749:   DM_DA *dd = (DM_DA *)da->data;

752:   if (refine_x) *refine_x = dd->refine_x;
753:   if (refine_y) *refine_y = dd->refine_y;
754:   if (refine_z) *refine_z = dd->refine_z;
755:   return 0;
756: }

758: /*@C
759:      DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.

761:     Logically Collective on da

763:   Input Parameters:
764: +    da - the DMDA object
765: -    f - the function that allocates the matrix for that specific DMDA

767:   Level: developer

769:    Notes:
770:     See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
771:        the diagonal and off-diagonal blocks of the matrix

773:    Not supported from Fortran

775: .seealso: `DMCreateMatrix()`, `DMDASetBlockFills()`
776: @*/
777: PetscErrorCode DMDASetGetMatrix(DM da, PetscErrorCode (*f)(DM, Mat *))
778: {
780:   da->ops->creatematrix = f;
781:   return 0;
782: }

784: /*@
785:    DMDAMapMatStencilToGlobal - Map a list of MatStencils on a grid to global indices.

787:    Not Collective

789:    Input Parameters:
790: +  da - the DMDA object
791: .  m - number of MatStencils
792: -  idxm - grid points (and component number when dof > 1)

794:    Output Parameter:
795: .  gidxm - global row indices

797:    Level: intermediate

799: .seealso: `MatStencil`
800: @*/
801: PetscErrorCode DMDAMapMatStencilToGlobal(DM da, PetscInt m, const MatStencil idxm[], PetscInt gidxm[])
802: {
803:   const DM_DA           *dd  = (const DM_DA *)da->data;
804:   const PetscInt        *dxm = (const PetscInt *)idxm;
805:   PetscInt               i, j, sdim, tmp, dim;
806:   PetscInt               dims[4], starts[4], dims2[3], starts2[3], dof = dd->w;
807:   ISLocalToGlobalMapping ltog;

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

811:   /* Code adapted from DMDAGetGhostCorners() */
812:   starts2[0] = dd->Xs / dof + dd->xo;
813:   starts2[1] = dd->Ys + dd->yo;
814:   starts2[2] = dd->Zs + dd->zo;
815:   dims2[0]   = (dd->Xe - dd->Xs) / dof;
816:   dims2[1]   = (dd->Ye - dd->Ys);
817:   dims2[2]   = (dd->Ze - dd->Zs);

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

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

843:   /* Map local indices to global indices */
844:   DMGetLocalToGlobalMapping(da, &ltog);
845:   ISLocalToGlobalMappingApply(ltog, m, gidxm, gidxm);
846:   return 0;
847: }

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

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

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

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

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

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

928: PetscErrorCode DMRefine_DA(DM da, MPI_Comm comm, DM *daref)
929: {
930:   PetscInt M, N, P, i, dim;
931:   Vec      coordsc, coordsf;
932:   DM       da2;
933:   DM_DA   *dd = (DM_DA *)da->data, *dd2;


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

995:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
996:   da2->ops->creatematrix = da->ops->creatematrix;
997:   /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
998:   da2->ops->getcoloring = da->ops->getcoloring;
999:   dd2->interptype       = dd->interptype;

1001:   /* copy fill information if given */
1002:   if (dd->dfill) {
1003:     PetscMalloc1(dd->dfill[dd->w] + dd->w + 1, &dd2->dfill);
1004:     PetscArraycpy(dd2->dfill, dd->dfill, dd->dfill[dd->w] + dd->w + 1);
1005:   }
1006:   if (dd->ofill) {
1007:     PetscMalloc1(dd->ofill[dd->w] + dd->w + 1, &dd2->ofill);
1008:     PetscArraycpy(dd2->ofill, dd->ofill, dd->ofill[dd->w] + dd->w + 1);
1009:   }
1010:   /* copy the refine information */
1011:   dd2->coarsen_x = dd2->refine_x = dd->refine_x;
1012:   dd2->coarsen_y = dd2->refine_y = dd->refine_y;
1013:   dd2->coarsen_z = dd2->refine_z = dd->refine_z;

1015:   if (dd->refine_z_hier) {
1016:     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];
1017:     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];
1018:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1019:     PetscMalloc1(dd2->refine_z_hier_n, &dd2->refine_z_hier);
1020:     PetscArraycpy(dd2->refine_z_hier, dd->refine_z_hier, dd2->refine_z_hier_n);
1021:   }
1022:   if (dd->refine_y_hier) {
1023:     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];
1024:     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];
1025:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1026:     PetscMalloc1(dd2->refine_y_hier_n, &dd2->refine_y_hier);
1027:     PetscArraycpy(dd2->refine_y_hier, dd->refine_y_hier, dd2->refine_y_hier_n);
1028:   }
1029:   if (dd->refine_x_hier) {
1030:     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];
1031:     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];
1032:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1033:     PetscMalloc1(dd2->refine_x_hier_n, &dd2->refine_x_hier);
1034:     PetscArraycpy(dd2->refine_x_hier, dd->refine_x_hier, dd2->refine_x_hier_n);
1035:   }

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

1040:   dd2->lf = dd->lf;
1041:   dd2->lj = dd->lj;

1043:   da2->leveldown = da->leveldown;
1044:   da2->levelup   = da->levelup + 1;

1046:   DMSetUp(da2);

1048:   /* interpolate coordinates if they are set on the coarse grid */
1049:   DMGetCoordinates(da, &coordsc);
1050:   if (coordsc) {
1051:     DM  cdaf, cdac;
1052:     Mat II;

1054:     DMGetCoordinateDM(da, &cdac);
1055:     DMGetCoordinateDM(da2, &cdaf);
1056:     /* force creation of the coordinate vector */
1057:     DMDASetUniformCoordinates(da2, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
1058:     DMGetCoordinates(da2, &coordsf);
1059:     DMCreateInterpolation(cdac, cdaf, &II, NULL);
1060:     MatInterpolate(II, coordsc, coordsf);
1061:     MatDestroy(&II);
1062:   }

1064:   for (i = 0; i < da->bs; i++) {
1065:     const char *fieldname;
1066:     DMDAGetFieldName(da, i, &fieldname);
1067:     DMDASetFieldName(da2, i, fieldname);
1068:   }

1070:   *daref = da2;
1071:   return 0;
1072: }

1074: PetscErrorCode DMCoarsen_DA(DM dmf, MPI_Comm comm, DM *dmc)
1075: {
1076:   PetscInt M, N, P, i, dim;
1077:   Vec      coordsc, coordsf;
1078:   DM       dmc2;
1079:   DM_DA   *dd = (DM_DA *)dmf->data, *dd2;


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

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

1147:   /* copy fill information if given */
1148:   if (dd->dfill) {
1149:     PetscMalloc1(dd->dfill[dd->w] + dd->w + 1, &dd2->dfill);
1150:     PetscArraycpy(dd2->dfill, dd->dfill, dd->dfill[dd->w] + dd->w + 1);
1151:   }
1152:   if (dd->ofill) {
1153:     PetscMalloc1(dd->ofill[dd->w] + dd->w + 1, &dd2->ofill);
1154:     PetscArraycpy(dd2->ofill, dd->ofill, dd->ofill[dd->w] + dd->w + 1);
1155:   }
1156:   /* copy the refine information */
1157:   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1158:   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1159:   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;

1161:   if (dd->refine_z_hier) {
1162:     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];
1163:     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];
1164:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1165:     PetscMalloc1(dd2->refine_z_hier_n, &dd2->refine_z_hier);
1166:     PetscArraycpy(dd2->refine_z_hier, dd->refine_z_hier, dd2->refine_z_hier_n);
1167:   }
1168:   if (dd->refine_y_hier) {
1169:     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];
1170:     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];
1171:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1172:     PetscMalloc1(dd2->refine_y_hier_n, &dd2->refine_y_hier);
1173:     PetscArraycpy(dd2->refine_y_hier, dd->refine_y_hier, dd2->refine_y_hier_n);
1174:   }
1175:   if (dd->refine_x_hier) {
1176:     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];
1177:     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];
1178:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1179:     PetscMalloc1(dd2->refine_x_hier_n, &dd2->refine_x_hier);
1180:     PetscArraycpy(dd2->refine_x_hier, dd->refine_x_hier, dd2->refine_x_hier_n);
1181:   }

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

1186:   dd2->lf = dd->lf;
1187:   dd2->lj = dd->lj;

1189:   dmc2->leveldown = dmf->leveldown + 1;
1190:   dmc2->levelup   = dmf->levelup;

1192:   DMSetUp(dmc2);

1194:   /* inject coordinates if they are set on the fine grid */
1195:   DMGetCoordinates(dmf, &coordsf);
1196:   if (coordsf) {
1197:     DM         cdaf, cdac;
1198:     Mat        inject;
1199:     VecScatter vscat;

1201:     DMGetCoordinateDM(dmf, &cdaf);
1202:     DMGetCoordinateDM(dmc2, &cdac);
1203:     /* force creation of the coordinate vector */
1204:     DMDASetUniformCoordinates(dmc2, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
1205:     DMGetCoordinates(dmc2, &coordsc);

1207:     DMCreateInjection(cdac, cdaf, &inject);
1208:     MatScatterGetVecScatter(inject, &vscat);
1209:     VecScatterBegin(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD);
1210:     VecScatterEnd(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD);
1211:     MatDestroy(&inject);
1212:   }

1214:   for (i = 0; i < dmf->bs; i++) {
1215:     const char *fieldname;
1216:     DMDAGetFieldName(dmf, i, &fieldname);
1217:     DMDASetFieldName(dmc2, i, fieldname);
1218:   }

1220:   *dmc = dmc2;
1221:   return 0;
1222: }

1224: PetscErrorCode DMRefineHierarchy_DA(DM da, PetscInt nlevels, DM daf[])
1225: {
1226:   PetscInt i, n, *refx, *refy, *refz;

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

1233:   /* Get refinement factors, defaults taken from the coarse DMDA */
1234:   PetscMalloc3(nlevels, &refx, nlevels, &refy, nlevels, &refz);
1235:   for (i = 0; i < nlevels; i++) DMDAGetRefinementFactor(da, &refx[i], &refy[i], &refz[i]);
1236:   n = nlevels;
1237:   PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_x", refx, &n, NULL);
1238:   n = nlevels;
1239:   PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_y", refy, &n, NULL);
1240:   n = nlevels;
1241:   PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_z", refz, &n, NULL);

1243:   DMDASetRefinementFactor(da, refx[0], refy[0], refz[0]);
1244:   DMRefine(da, PetscObjectComm((PetscObject)da), &daf[0]);
1245:   for (i = 1; i < nlevels; i++) {
1246:     DMDASetRefinementFactor(daf[i - 1], refx[i], refy[i], refz[i]);
1247:     DMRefine(daf[i - 1], PetscObjectComm((PetscObject)da), &daf[i]);
1248:   }
1249:   PetscFree3(refx, refy, refz);
1250:   return 0;
1251: }

1253: PetscErrorCode DMCoarsenHierarchy_DA(DM da, PetscInt nlevels, DM dac[])
1254: {
1255:   PetscInt i;

1259:   if (nlevels == 0) return 0;
1261:   DMCoarsen(da, PetscObjectComm((PetscObject)da), &dac[0]);
1262:   for (i = 1; i < nlevels; i++) DMCoarsen(dac[i - 1], PetscObjectComm((PetscObject)da), &dac[i]);
1263:   return 0;
1264: }

1266: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm, PetscInt n, PetscReal *nodes)
1267: {
1268:   PetscInt     i, j, xs, xn, q;
1269:   PetscScalar *xx;
1270:   PetscReal    h;
1271:   Vec          x;
1272:   DM_DA       *da = (DM_DA *)dm->data;

1274:   if (da->bx != DM_BOUNDARY_PERIODIC) {
1275:     DMDAGetInfo(dm, NULL, &q, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
1276:     q = (q - 1) / (n - 1); /* number of spectral elements */
1277:     h = 2.0 / q;
1278:     DMDAGetCorners(dm, &xs, NULL, NULL, &xn, NULL, NULL);
1279:     xs = xs / (n - 1);
1280:     xn = xn / (n - 1);
1281:     DMDASetUniformCoordinates(dm, -1., 1., 0., 0., 0., 0.);
1282:     DMGetCoordinates(dm, &x);
1283:     DMDAVecGetArray(dm, x, &xx);

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

1297: /*@

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

1301:    Collective on da

1303:    Input Parameters:
1304: +   da - the DMDA object
1305: -   n - the number of GLL nodes
1306: -   nodes - the GLL nodes

1308:    Notes:
1309:     the parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1310:           on each process much be divisible by the number of GLL elements needed per process. This depends on whether the DM is
1311:           periodic or not.

1313:    Level: advanced

1315: .seealso: `DMDACreate()`, `PetscDTGaussLobattoLegendreQuadrature()`, `DMGetCoordinates()`
1316: @*/
1317: PetscErrorCode DMDASetGLLCoordinates(DM da, PetscInt n, PetscReal *nodes)
1318: {
1319:   if (da->dim == 1) {
1320:     DMDASetGLLCoordinates_1d(da, n, nodes);
1321:   } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not yet implemented for 2 or 3d");
1322:   return 0;
1323: }

1325: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM da1, DM dm2, PetscBool *compatible, PetscBool *set)
1326: {
1327:   DM_DA    *dd1 = (DM_DA *)da1->data, *dd2;
1328:   DM        da2;
1329:   DMType    dmtype2;
1330:   PetscBool isda, compatibleLocal;
1331:   PetscInt  i;

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