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, <og);
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: }