Actual source code: da.c
1: #include <petsc/private/dmdaimpl.h>
3: /*@
4: DMDASetSizes - Sets the number of grid points in the three dimensional directions
6: Logically Collective
8: Input Parameters:
9: + da - the `DMDA`
10: . M - the global X size
11: . N - the global Y size
12: - P - the global Z size
14: Level: intermediate
16: Developer Note:
17: Since the dimension may not yet have been set the code cannot error check for non-positive Y and Z number of grid points
19: .seealso: [](sec_struct), `DM`, `DMDA`, `PetscSplitOwnership()`
20: @*/
21: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
22: {
23: DM_DA *dd = (DM_DA *)da->data;
25: PetscFunctionBegin;
30: PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
31: PetscCheck(M >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Number of grid points in X direction must be positive");
32: PetscCheck(N >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Number of grid points in Y direction must be positive");
33: PetscCheck(P >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Number of grid points in Z direction must be positive");
35: dd->M = M;
36: dd->N = N;
37: dd->P = P;
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: /*@
42: DMDASetNumProcs - Sets the number of processes in each dimension
44: Logically Collective
46: Input Parameters:
47: + da - the `DMDA`
48: . m - the number of X processes (or `PETSC_DECIDE`)
49: . n - the number of Y processes (or `PETSC_DECIDE`)
50: - p - the number of Z processes (or `PETSC_DECIDE`)
52: Level: intermediate
54: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDASetSizes()`, `PetscSplitOwnership()`
55: @*/
56: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
57: {
58: DM_DA *dd = (DM_DA *)da->data;
60: PetscFunctionBegin;
65: PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
66: dd->m = m;
67: dd->n = n;
68: dd->p = p;
69: if (da->dim == 2) {
70: PetscMPIInt size;
71: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
72: if ((dd->m > 0) && (dd->n < 0)) {
73: dd->n = size / dd->m;
74: PetscCheck(dd->n * dd->m == size, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " processes in X direction not divisible into comm size %d", m, size);
75: }
76: if ((dd->n > 0) && (dd->m < 0)) {
77: dd->m = size / dd->n;
78: PetscCheck(dd->n * dd->m == size, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " processes in Y direction not divisible into comm size %d", n, size);
79: }
80: }
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: /*@
85: DMDAGetBoundaryType - Gets the type of ghost nodes on domain boundaries.
87: Not Collective
89: Input Parameter:
90: . da - The `DMDA`
92: Output Parameters:
93: + bx - x boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
94: . by - y boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
95: - bz - z boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
97: Level: intermediate
99: .seealso: [](sec_struct), `DMDASetBoundaryType()`, `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMBoundaryType`, `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
100: @*/
101: PetscErrorCode DMDAGetBoundaryType(DM da, PeOp DMBoundaryType *bx, PeOp DMBoundaryType *by, PeOp DMBoundaryType *bz)
102: {
103: DM_DA *dd = (DM_DA *)da->data;
105: PetscFunctionBegin;
107: if (bx) PetscAssertPointer(bx, 2);
108: if (by) PetscAssertPointer(by, 3);
109: if (bz) PetscAssertPointer(bz, 4);
110: if (bx) *bx = dd->bx;
111: if (by) *by = dd->by;
112: if (bz) *bz = dd->bz;
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*@
117: DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries for a `DMDA` object.
119: Not Collective
121: Input Parameters:
122: + da - The `DMDA`
123: . bx - x boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
124: . by - y boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
125: - bz - z boundary type, one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
127: Level: intermediate
129: Note:
130: The default is `DM_BOUNDARY_NONE`
132: .seealso: [](sec_struct), `DMDAGetBoundaryType()`, `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMBoundaryType`, `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`
133: @*/
134: PetscErrorCode DMDASetBoundaryType(DM da, DMBoundaryType bx, DMBoundaryType by, DMBoundaryType bz)
135: {
136: DM_DA *dd = (DM_DA *)da->data;
138: PetscFunctionBegin;
143: PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
144: dd->bx = bx;
145: dd->by = by;
146: dd->bz = bz;
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: /*@
151: DMDASetDof - Sets the number of degrees of freedom per vertex
153: Not Collective
155: Input Parameters:
156: + da - The `DMDA`
157: - dof - Number of degrees of freedom per vertex
159: Level: intermediate
161: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetDof()`, `DMDACreate()`, `DMDestroy()`
162: @*/
163: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
164: {
165: DM_DA *dd = (DM_DA *)da->data;
167: PetscFunctionBegin;
170: PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
171: dd->w = dof;
172: da->bs = dof;
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*@
177: DMDAGetDof - Gets the number of degrees of freedom per vertex
179: Not Collective
181: Input Parameter:
182: . da - The `DMDA`
184: Output Parameter:
185: . dof - Number of degrees of freedom per vertex
187: Level: intermediate
189: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDASetDof()`, `DMDACreate()`, `DMDestroy()`
190: @*/
191: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
192: {
193: DM_DA *dd = (DM_DA *)da->data;
195: PetscFunctionBegin;
197: PetscAssertPointer(dof, 2);
198: *dof = dd->w;
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: /*@
203: DMDAGetOverlap - Gets the size of the per-processor overlap.
205: Not Collective
207: Input Parameter:
208: . da - The `DMDA`
210: Output Parameters:
211: + x - Overlap in the x direction
212: . y - Overlap in the y direction
213: - z - Overlap in the z direction
215: Level: intermediate
217: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDASetOverlap()`
218: @*/
219: PetscErrorCode DMDAGetOverlap(DM da, PeOp PetscInt *x, PeOp PetscInt *y, PeOp PetscInt *z)
220: {
221: DM_DA *dd = (DM_DA *)da->data;
223: PetscFunctionBegin;
225: if (x) *x = dd->xol;
226: if (y) *y = dd->yol;
227: if (z) *z = dd->zol;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: DMDASetOverlap - Sets the size of the per-processor overlap.
234: Not Collective
236: Input Parameters:
237: + da - The `DMDA`
238: . x - Overlap in the x direction
239: . y - Overlap in the y direction
240: - z - Overlap in the z direction
242: Level: intermediate
244: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDAGetOverlap()`
245: @*/
246: PetscErrorCode DMDASetOverlap(DM da, PetscInt x, PetscInt y, PetscInt z)
247: {
248: DM_DA *dd = (DM_DA *)da->data;
250: PetscFunctionBegin;
255: dd->xol = x;
256: dd->yol = y;
257: dd->zol = z;
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /*@
262: DMDAGetNumLocalSubDomains - Gets the number of local subdomains that would be created upon decomposition.
264: Not Collective
266: Input Parameter:
267: . da - The `DMDA`
269: Output Parameter:
270: . Nsub - Number of local subdomains created upon decomposition
272: Level: intermediate
274: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDASetNumLocalSubDomains()`
275: @*/
276: PetscErrorCode DMDAGetNumLocalSubDomains(DM da, PetscInt *Nsub)
277: {
278: DM_DA *dd = (DM_DA *)da->data;
280: PetscFunctionBegin;
282: if (Nsub) *Nsub = dd->Nsub;
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@
287: DMDASetNumLocalSubDomains - Sets the number of local subdomains to create when decomposing with `DMCreateDomainDecomposition()`
289: Not Collective
291: Input Parameters:
292: + da - The `DMDA`
293: - Nsub - The number of local subdomains requested
295: Level: intermediate
297: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateDomainDecomposition()`, `DMDAGetNumLocalSubDomains()`
298: @*/
299: PetscErrorCode DMDASetNumLocalSubDomains(DM da, PetscInt Nsub)
300: {
301: DM_DA *dd = (DM_DA *)da->data;
303: PetscFunctionBegin;
306: dd->Nsub = Nsub;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*@
311: DMDASetOffset - Sets the index offset of the `DMDA`.
313: Collective
315: Input Parameters:
316: + da - The `DMDA`
317: . xo - The offset in the x direction
318: . yo - The offset in the y direction
319: . zo - The offset in the z direction
320: . Mo - The problem offset in the x direction
321: . No - The problem offset in the y direction
322: - Po - The problem offset in the z direction
324: Level: developer
326: Note:
327: This is used primarily to overlap a computation on a local `DMDA` with that on a global `DMDA` without
328: changing boundary conditions or subdomain features that depend upon the global offsets.
330: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
331: @*/
332: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
333: {
334: DM_DA *dd = (DM_DA *)da->data;
336: PetscFunctionBegin;
344: dd->xo = xo;
345: dd->yo = yo;
346: dd->zo = zo;
347: dd->Mo = Mo;
348: dd->No = No;
349: dd->Po = Po;
351: if (da->coordinates[0].dm) PetscCall(DMDASetOffset(da->coordinates[0].dm, xo, yo, zo, Mo, No, Po));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@
356: DMDAGetOffset - Gets the index offset of the `DMDA`.
358: Not Collective
360: Input Parameter:
361: . da - The `DMDA`
363: Output Parameters:
364: + xo - The offset in the x direction
365: . yo - The offset in the y direction
366: . zo - The offset in the z direction
367: . Mo - The global size in the x direction
368: . No - The global size in the y direction
369: - Po - The global size in the z direction
371: Level: developer
373: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDASetOffset()`, `DMDAVecGetArray()`
374: @*/
375: PetscErrorCode DMDAGetOffset(DM da, PeOp PetscInt *xo, PeOp PetscInt *yo, PeOp PetscInt *zo, PeOp PetscInt *Mo, PeOp PetscInt *No, PeOp PetscInt *Po)
376: {
377: DM_DA *dd = (DM_DA *)da->data;
379: PetscFunctionBegin;
381: if (xo) *xo = dd->xo;
382: if (yo) *yo = dd->yo;
383: if (zo) *zo = dd->zo;
384: if (Mo) *Mo = dd->Mo;
385: if (No) *No = dd->No;
386: if (Po) *Po = dd->Po;
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@
391: DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain `DMDA`.
393: Not Collective
395: Input Parameter:
396: . da - The `DMDA`
398: Output Parameters:
399: + xs - The start of the region in x
400: . ys - The start of the region in y
401: . zs - The start of the region in z
402: . xm - The size of the region in x
403: . ym - The size of the region in y
404: - zm - The size of the region in z
406: Level: intermediate
408: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
409: @*/
410: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PeOp PetscInt *xs, PeOp PetscInt *ys, PeOp PetscInt *zs, PeOp PetscInt *xm, PeOp PetscInt *ym, PeOp PetscInt *zm)
411: {
412: DM_DA *dd = (DM_DA *)da->data;
414: PetscFunctionBegin;
416: if (xs) *xs = dd->nonxs;
417: if (ys) *ys = dd->nonys;
418: if (zs) *zs = dd->nonzs;
419: if (xm) *xm = dd->nonxm;
420: if (ym) *ym = dd->nonym;
421: if (zm) *zm = dd->nonzm;
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*@
426: DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain `DMDA`.
428: Collective
430: Input Parameters:
431: + da - The `DMDA`
432: . xs - The start of the region in x
433: . ys - The start of the region in y
434: . zs - The start of the region in z
435: . xm - The size of the region in x
436: . ym - The size of the region in y
437: - zm - The size of the region in z
439: Level: intermediate
441: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetOffset()`, `DMDAVecGetArray()`
442: @*/
443: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
444: {
445: DM_DA *dd = (DM_DA *)da->data;
447: PetscFunctionBegin;
455: dd->nonxs = xs;
456: dd->nonys = ys;
457: dd->nonzs = zs;
458: dd->nonxm = xm;
459: dd->nonym = ym;
460: dd->nonzm = zm;
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@
465: DMDASetStencilType - Sets the type of the communication stencil
467: Logically Collective
469: Input Parameters:
470: + da - The `DMDA`
471: - stype - The stencil type, use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.
473: Level: intermediate
475: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMDAStencilType`, `DMDA_STENCIL_BOX`, `DMDA_STENCIL_STAR.`
476: @*/
477: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
478: {
479: DM_DA *dd = (DM_DA *)da->data;
481: PetscFunctionBegin;
484: PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
485: dd->stencil_type = stype;
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: /*@
490: DMDAGetStencilType - Gets the type of the communication stencil
492: Not Collective
494: Input Parameter:
495: . da - The `DMDA`
497: Output Parameter:
498: . stype - The stencil type, use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`.
500: Level: intermediate
502: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMDAStencilType`, `DMDA_STENCIL_BOX`, `DMDA_STENCIL_STAR.`
503: @*/
504: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
505: {
506: DM_DA *dd = (DM_DA *)da->data;
508: PetscFunctionBegin;
510: PetscAssertPointer(stype, 2);
511: *stype = dd->stencil_type;
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: /*@
516: DMDASetStencilWidth - Sets the width of the communication stencil
518: Logically Collective
520: Input Parameters:
521: + da - The `DMDA`
522: - width - The stencil width
524: Level: intermediate
526: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMDAStencilType`, `DMDA_STENCIL_BOX`, `DMDA_STENCIL_STAR.`
527: @*/
528: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
529: {
530: DM_DA *dd = (DM_DA *)da->data;
532: PetscFunctionBegin;
535: PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
536: dd->s = width;
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: /*@
541: DMDAGetStencilWidth - Gets the width of the communication stencil
543: Not Collective
545: Input Parameter:
546: . da - The `DMDA`
548: Output Parameter:
549: . width - The stencil width
551: Level: intermediate
553: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`, `DMDAStencilType`, `DMDA_STENCIL_BOX`, `DMDA_STENCIL_STAR.`
554: @*/
555: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
556: {
557: DM_DA *dd = (DM_DA *)da->data;
559: PetscFunctionBegin;
561: PetscAssertPointer(width, 2);
562: *width = dd->s;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da, PetscInt M, PetscInt m, const PetscInt lx[])
567: {
568: PetscInt i, sum;
570: PetscFunctionBegin;
571: PetscCheck(M >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Global dimension not set");
572: for (i = sum = 0; i < m; i++) sum += lx[i];
573: PetscCheck(sum == M, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_INCOMP, "Ownership ranges sum to %" PetscInt_FMT " but global dimension is %" PetscInt_FMT, sum, M);
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
580: Logically Collective
582: Input Parameters:
583: + da - The `DMDA`
584: . lx - array containing number of nodes in the X direction on each process, or `NULL`. If non-null, must be of length da->m
585: . ly - array containing number of nodes in the Y direction on each process, or `NULL`. If non-null, must be of length da->n
586: - lz - array containing number of nodes in the Z direction on each process, or `NULL`. If non-null, must be of length da->p.
588: Level: intermediate
590: Note:
591: These numbers are NOT multiplied by the number of dof per node.
593: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `DMDestroy()`
594: @*/
595: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
596: {
597: DM_DA *dd = (DM_DA *)da->data;
599: PetscFunctionBegin;
601: PetscCheck(!da->setupcalled, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "This function must be called before DMSetUp()");
602: if (lx) {
603: PetscCheck(dd->m >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Cannot set ownership ranges before setting number of processes");
604: PetscCall(DMDACheckOwnershipRanges_Private(da, dd->M, dd->m, lx));
605: if (!dd->lx) PetscCall(PetscMalloc1(dd->m, &dd->lx));
606: PetscCall(PetscArraycpy(dd->lx, lx, dd->m));
607: }
608: if (ly) {
609: PetscCheck(dd->n >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Cannot set ownership ranges before setting number of processes");
610: PetscCall(DMDACheckOwnershipRanges_Private(da, dd->N, dd->n, ly));
611: if (!dd->ly) PetscCall(PetscMalloc1(dd->n, &dd->ly));
612: PetscCall(PetscArraycpy(dd->ly, ly, dd->n));
613: }
614: if (lz) {
615: PetscCheck(dd->p >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONGSTATE, "Cannot set ownership ranges before setting number of processes");
616: PetscCall(DMDACheckOwnershipRanges_Private(da, dd->P, dd->p, lz));
617: if (!dd->lz) PetscCall(PetscMalloc1(dd->p, &dd->lz));
618: PetscCall(PetscArraycpy(dd->lz, lz, dd->p));
619: }
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: /*@
624: DMDASetInterpolationType - Sets the type of interpolation that will be
625: returned by `DMCreateInterpolation()`
627: Logically Collective
629: Input Parameters:
630: + da - initial distributed array
631: - ctype - `DMDA_Q1` and `DMDA_Q0` are currently the only supported forms
633: Level: intermediate
635: Note:
636: You should call this on the coarser of the two `DMDA` you pass to `DMCreateInterpolation()`
638: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMDestroy()`, `DMDAInterpolationType`,
639: `DMDA_Q1`, `DMDA_Q0`
640: @*/
641: PetscErrorCode DMDASetInterpolationType(DM da, DMDAInterpolationType ctype)
642: {
643: DM_DA *dd = (DM_DA *)da->data;
645: PetscFunctionBegin;
648: dd->interptype = ctype;
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: /*@
653: DMDAGetInterpolationType - Gets the type of interpolation that will be
654: used by `DMCreateInterpolation()`
656: Not Collective
658: Input Parameter:
659: . da - distributed array
661: Output Parameter:
662: . ctype - interpolation type (`DMDA_Q1` and `DMDA_Q0` are currently the only supported forms)
664: Level: intermediate
666: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAInterpolationType`, `DMDASetInterpolationType()`, `DMCreateInterpolation()`,
667: `DMDA_Q1`, `DMDA_Q0`
668: @*/
669: PetscErrorCode DMDAGetInterpolationType(DM da, DMDAInterpolationType *ctype)
670: {
671: DM_DA *dd = (DM_DA *)da->data;
673: PetscFunctionBegin;
675: PetscAssertPointer(ctype, 2);
676: *ctype = dd->interptype;
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: /*@C
681: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
682: processes neighbors.
684: Not Collective
686: Input Parameter:
687: . da - the `DMDA` object
689: Output Parameter:
690: . ranks - the neighbors ranks, stored with the x index increasing most rapidly. The process itself is in the list
692: Level: intermediate
694: Notes:
695: In 2d the `ranks` is of length 9, in 3d of length 27
697: Not supported in 1d
699: Do not free the array, it is freed when the `DMDA` is destroyed.
701: Fortran Note:
702: Use `DMDARestoreNeighbors()` to return the array when no longer needed
704: .seealso: [](sec_struct), `DMDA`, `DM`
705: @*/
706: PetscErrorCode DMDAGetNeighbors(DM da, const PetscMPIInt *ranks[])
707: {
708: DM_DA *dd = (DM_DA *)da->data;
710: PetscFunctionBegin;
712: *ranks = dd->neighbors;
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: /*@C
717: DMDAGetOwnershipRanges - Gets the number of indices in the x, y and z direction that are owned by each process in that direction
719: Not Collective
721: Input Parameter:
722: . da - the `DMDA` object
724: Output Parameters:
725: + lx - ownership along x direction (optional), its length is `m` the number of processes in the x-direction
726: . ly - ownership along y direction (optional), its length is `n` the number of processes in the y-direction
727: - lz - ownership along z direction (optional), its length is `p` the number of processes in the z-direction
729: Level: intermediate
731: Notes:
732: These correspond to the optional final arguments passed to `DMDACreate()`, `DMDACreate2d()`, `DMDACreate3d()`
734: You should not free these arrays, nor change the values in them. They will only have valid values while the
735: `DMDA` they came from still exists (has not been destroyed).
737: These numbers are NOT multiplied by the number of dof per node.
739: The meaning of these is different than that returned by `VecGetOwnerShipRanges()`
741: Fortran Notes:
742: Pass `PETSC_NULL_INT_POINTER` for any array not needed.
744: Use `DMDARestoreOwershipRange()` to return the arrays when no longer needed
746: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetCorners()`, `DMDAGetGhostCorners()`, `DMDACreate()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`, `VecGetOwnershipRanges()`
747: @*/
748: PetscErrorCode DMDAGetOwnershipRanges(DM da, PeOp const PetscInt *lx[], PeOp const PetscInt *ly[], PeOp const PetscInt *lz[])
749: {
750: DM_DA *dd = (DM_DA *)da->data;
752: PetscFunctionBegin;
754: if (lx) *lx = dd->lx;
755: if (ly) *ly = dd->ly;
756: if (lz) *lz = dd->lz;
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: /*@
761: DMDASetRefinementFactor - Set the ratios that the `DMDA` grid is refined
763: Logically Collective
765: Input Parameters:
766: + da - the `DMDA` object
767: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
768: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
769: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
771: Options Database Keys:
772: + -da_refine_x refine_x - refinement ratio in x direction
773: . -da_refine_y rafine_y - refinement ratio in y direction
774: . -da_refine_z refine_z - refinement ratio in z direction
775: - -da_refine <n> - refine the `DMDA` object n times when it is created.
777: Level: intermediate
779: Note:
780: Pass `PETSC_IGNORE` to leave a value unchanged
782: .seealso: [](sec_struct), `DM`, `DMDA`, `DMRefine()`, `DMDAGetRefinementFactor()`
783: @*/
784: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y, PetscInt refine_z)
785: {
786: DM_DA *dd = (DM_DA *)da->data;
788: PetscFunctionBegin;
794: if (refine_x > 0) dd->refine_x = refine_x;
795: if (refine_y > 0) dd->refine_y = refine_y;
796: if (refine_z > 0) dd->refine_z = refine_z;
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }
800: /*@
801: DMDAGetRefinementFactor - Gets the ratios that the `DMDA` grid is refined
803: Not Collective
805: Input Parameter:
806: . da - the `DMDA` object
808: Output Parameters:
809: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
810: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
811: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
813: Level: intermediate
815: Note:
816: Pass `NULL` for values you do not need
818: .seealso: [](sec_struct), `DM`, `DMDA`, `DMRefine()`, `DMDASetRefinementFactor()`
819: @*/
820: PetscErrorCode DMDAGetRefinementFactor(DM da, PeOp PetscInt *refine_x, PeOp PetscInt *refine_y, PeOp PetscInt *refine_z)
821: {
822: DM_DA *dd = (DM_DA *)da->data;
824: PetscFunctionBegin;
826: if (refine_x) *refine_x = dd->refine_x;
827: if (refine_y) *refine_y = dd->refine_y;
828: if (refine_z) *refine_z = dd->refine_z;
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: /*@C
833: DMDASetGetMatrix - Sets the routine used by the `DMDA` to allocate a matrix.
835: Logically Collective; No Fortran Support
837: Input Parameters:
838: + da - the `DMDA` object
839: - f - the function that allocates the matrix for that specific `DMDA`
841: Calling sequence of `f`:
842: + da - the `DMDA` object
843: - A - the created matrix
845: Level: developer
847: Notes:
848: If the function is not provided a default function is used that uses the `DMDAStencilType`, `DMBoundaryType`, and value of `DMDASetStencilWidth()`
849: to construct the matrix.
851: See `DMDASetBlockFills()` that provides a simple way to provide the nonzero structure for
852: the diagonal and off-diagonal blocks of the matrix without providing a custom function
854: Developer Note:
855: This should be called `DMDASetCreateMatrix()`
857: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateMatrix()`, `DMDASetBlockFills()`
858: @*/
859: PetscErrorCode DMDASetGetMatrix(DM da, PetscErrorCode (*f)(DM da, Mat *A))
860: {
861: PetscFunctionBegin;
863: da->ops->creatematrix = f;
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: DMDAMapMatStencilToGlobal - Map a list of `MatStencil` on a grid to global indices.
870: Not Collective
872: Input Parameters:
873: + da - the `DMDA` object
874: . m - number of `MatStencil` to map
875: - idxm - grid points (and component number when dof > 1)
877: Output Parameter:
878: . gidxm - global row indices
880: Level: intermediate
882: .seealso: [](sec_struct), `DM`, `DMDA`, `MatStencil`
883: @*/
884: PetscErrorCode DMDAMapMatStencilToGlobal(DM da, PetscInt m, const MatStencil idxm[], PetscInt gidxm[])
885: {
886: const DM_DA *dd = (const DM_DA *)da->data;
887: const PetscInt *dxm = (const PetscInt *)idxm;
888: PetscInt i, j, sdim, tmp, dim;
889: PetscInt dims[4], starts[4], dims2[3], starts2[3], dof = dd->w;
890: ISLocalToGlobalMapping ltog;
892: PetscFunctionBegin;
893: if (m <= 0) PetscFunctionReturn(PETSC_SUCCESS);
895: /* Code adapted from DMDAGetGhostCorners() */
896: starts2[0] = dd->Xs / dof + dd->xo;
897: starts2[1] = dd->Ys + dd->yo;
898: starts2[2] = dd->Zs + dd->zo;
899: dims2[0] = (dd->Xe - dd->Xs) / dof;
900: dims2[1] = (dd->Ye - dd->Ys);
901: dims2[2] = (dd->Ze - dd->Zs);
903: /* As if we do MatSetStencil() to get dims[]/starts[] of mat->stencil */
904: dim = da->dim; /* DA dim: 1 to 3 */
905: sdim = dim + (dof > 1 ? 1 : 0); /* Dimensions in MatStencil's (k,j,i,c) view */
906: for (i = 0; i < dim; i++) { /* Reverse the order and also skip the unused dimensions */
907: dims[i] = dims2[dim - i - 1]; /* ex. dims/starts[] are in order of {i} for 1D, {j,i} for 2D and {k,j,i} for 3D */
908: starts[i] = starts2[dim - i - 1];
909: }
910: starts[dim] = 0; /* Append the extra dim for dof (won't be used below if dof=1) */
911: dims[dim] = dof;
913: /* Map stencils to local indices (code adapted from MatSetValuesStencil()) */
914: for (i = 0; i < m; i++) {
915: dxm += 3 - dim; /* Input is {k,j,i,c}; move the pointer to the first used index, e.g., j in 2D */
916: tmp = 0;
917: for (j = 0; j < sdim; j++) { /* Iter over, ex. j,i or j,i,c in 2D */
918: if (tmp < 0 || dxm[j] < starts[j] || dxm[j] >= (starts[j] + dims[j])) tmp = -1; /* Beyond the ghost region, therefore ignored with negative indices */
919: else tmp = tmp * dims[j] + (dxm[j] - starts[j]);
920: }
921: gidxm[i] = tmp;
922: /* Move to the next MatStencil point */
923: if (dof > 1) dxm += sdim; /* c is already counted in sdim */
924: else dxm += sdim + 1; /* skip the unused c */
925: }
927: /* Map local indices to global indices */
928: PetscCall(DMGetLocalToGlobalMapping(da, <og));
929: PetscCall(ISLocalToGlobalMappingApply(ltog, m, gidxm, gidxm));
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: /*
934: Creates "balanced" ownership ranges after refinement, constrained by the need for the
935: fine grid boundaries to fall within one stencil width of the coarse partition.
937: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
938: */
939: static PetscErrorCode DMDARefineOwnershipRanges(DM da, PetscBool periodic, PetscInt stencil_width, PetscInt ratio, PetscInt m, const PetscInt lc[], PetscInt lf[])
940: {
941: PetscInt i, totalc = 0, remaining, startc = 0, startf = 0;
943: PetscFunctionBegin;
944: PetscCheck(ratio >= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Requested refinement ratio %" PetscInt_FMT " must be at least 1", ratio);
945: if (ratio == 1) {
946: PetscCall(PetscArraycpy(lf, lc, m));
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
949: for (i = 0; i < m; i++) totalc += lc[i];
950: remaining = (!periodic) + ratio * (totalc - (!periodic));
951: for (i = 0; i < m; i++) {
952: PetscInt want = remaining / (m - i) + !!(remaining % (m - i));
953: if (i == m - 1) lf[i] = want;
954: else {
955: const PetscInt nextc = startc + lc[i];
956: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
957: * coarse stencil width of the first coarse node in the next subdomain. */
958: while ((startf + want) / ratio < nextc - stencil_width) want++;
959: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
960: * coarse stencil width of the last coarse node in the current subdomain. */
961: while ((startf + want - 1 + ratio - 1) / ratio > nextc - 1 + stencil_width) want--;
962: /* Make sure all constraints are satisfied */
963: if (want < 0 || want > remaining || ((startf + want) / ratio < nextc - stencil_width) || ((startf + want - 1 + ratio - 1) / ratio > nextc - 1 + stencil_width))
964: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Could not find a compatible refined ownership range");
965: }
966: lf[i] = want;
967: startc += lc[i];
968: startf += lf[i];
969: remaining -= lf[i];
970: }
971: PetscFunctionReturn(PETSC_SUCCESS);
972: }
974: /*
975: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
976: fine grid boundaries to fall within one stencil width of the coarse partition.
978: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
979: */
980: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da, PetscBool periodic, PetscInt stencil_width, PetscInt ratio, PetscInt m, const PetscInt lf[], PetscInt lc[])
981: {
982: PetscInt i, totalf, remaining, startc, startf;
984: PetscFunctionBegin;
985: PetscCheck(ratio >= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Requested refinement ratio %" PetscInt_FMT " must be at least 1", ratio);
986: if (ratio == 1) {
987: PetscCall(PetscArraycpy(lc, lf, m));
988: PetscFunctionReturn(PETSC_SUCCESS);
989: }
990: for (i = 0, totalf = 0; i < m; i++) totalf += lf[i];
991: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
992: for (i = 0, startc = 0, startf = 0; i < m; i++) {
993: PetscInt want = remaining / (m - i) + !!(remaining % (m - i));
994: if (i == m - 1) lc[i] = want;
995: else {
996: const PetscInt nextf = startf + lf[i];
997: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
998: * node is within one stencil width. */
999: while (nextf / ratio < startc + want - stencil_width) want--;
1000: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
1001: * fine node is within one stencil width. */
1002: while ((nextf - 1 + ratio - 1) / ratio > startc + want - 1 + stencil_width) want++;
1003: if (want < 0 || want > remaining || (nextf / ratio < startc + want - stencil_width) || ((nextf - 1 + ratio - 1) / ratio > startc + want - 1 + stencil_width))
1004: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_SIZ, "Could not find a compatible coarsened ownership range");
1005: }
1006: lc[i] = want;
1007: startc += lc[i];
1008: startf += lf[i];
1009: remaining -= lc[i];
1010: }
1011: PetscFunctionReturn(PETSC_SUCCESS);
1012: }
1014: PetscErrorCode DMRefine_DA(DM da, MPI_Comm comm, DM *daref)
1015: {
1016: PetscInt M, N, P, i, dim;
1017: Vec coordsc, coordsf;
1018: DM da2;
1019: DM_DA *dd = (DM_DA *)da->data, *dd2;
1021: PetscFunctionBegin;
1023: PetscAssertPointer(daref, 3);
1025: PetscCall(DMGetDimension(da, &dim));
1026: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1027: M = dd->refine_x * dd->M;
1028: } else {
1029: M = 1 + dd->refine_x * (dd->M - 1);
1030: }
1031: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1032: if (dim > 1) {
1033: N = dd->refine_y * dd->N;
1034: } else {
1035: N = 1;
1036: }
1037: } else {
1038: N = 1 + dd->refine_y * (dd->N - 1);
1039: }
1040: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1041: if (dim > 2) {
1042: P = dd->refine_z * dd->P;
1043: } else {
1044: P = 1;
1045: }
1046: } else {
1047: P = 1 + dd->refine_z * (dd->P - 1);
1048: }
1049: PetscCall(DMDACreate(PetscObjectComm((PetscObject)da), &da2));
1050: PetscCall(DMSetOptionsPrefix(da2, ((PetscObject)da)->prefix));
1051: PetscCall(DMSetDimension(da2, dim));
1052: PetscCall(DMDASetSizes(da2, M, N, P));
1053: PetscCall(DMDASetNumProcs(da2, dd->m, dd->n, dd->p));
1054: PetscCall(DMDASetBoundaryType(da2, dd->bx, dd->by, dd->bz));
1055: PetscCall(DMDASetDof(da2, dd->w));
1056: PetscCall(DMDASetStencilType(da2, dd->stencil_type));
1057: PetscCall(DMDASetStencilWidth(da2, dd->s));
1058: if (dim == 3) {
1059: PetscInt *lx, *ly, *lz;
1060: PetscCall(PetscMalloc3(dd->m, &lx, dd->n, &ly, dd->p, &lz));
1061: PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx));
1062: PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_y, dd->n, dd->ly, ly));
1063: PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_z, dd->p, dd->lz, lz));
1064: PetscCall(DMDASetOwnershipRanges(da2, lx, ly, lz));
1065: PetscCall(PetscFree3(lx, ly, lz));
1066: } else if (dim == 2) {
1067: PetscInt *lx, *ly;
1068: PetscCall(PetscMalloc2(dd->m, &lx, dd->n, &ly));
1069: PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx));
1070: PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_y, dd->n, dd->ly, ly));
1071: PetscCall(DMDASetOwnershipRanges(da2, lx, ly, NULL));
1072: PetscCall(PetscFree2(lx, ly));
1073: } else if (dim == 1) {
1074: PetscInt *lx;
1075: PetscCall(PetscMalloc1(dd->m, &lx));
1076: PetscCall(DMDARefineOwnershipRanges(da, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->refine_x, dd->m, dd->lx, lx));
1077: PetscCall(DMDASetOwnershipRanges(da2, lx, NULL, NULL));
1078: PetscCall(PetscFree(lx));
1079: }
1080: dd2 = (DM_DA *)da2->data;
1082: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
1083: da2->ops->creatematrix = da->ops->creatematrix;
1084: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
1085: da2->ops->getcoloring = da->ops->getcoloring;
1086: dd2->interptype = dd->interptype;
1088: /* copy fill information if given */
1089: if (dd->dfill) {
1090: PetscCall(PetscMalloc1(dd->dfill[dd->w] + dd->w + 1, &dd2->dfill));
1091: PetscCall(PetscArraycpy(dd2->dfill, dd->dfill, dd->dfill[dd->w] + dd->w + 1));
1092: }
1093: if (dd->ofill) {
1094: PetscCall(PetscMalloc1(dd->ofill[dd->w] + dd->w + 1, &dd2->ofill));
1095: PetscCall(PetscArraycpy(dd2->ofill, dd->ofill, dd->ofill[dd->w] + dd->w + 1));
1096: }
1097: /* copy the refine information */
1098: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
1099: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
1100: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
1102: if (dd->refine_z_hier) {
1103: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
1104: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
1105: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1106: PetscCall(PetscMalloc1(dd2->refine_z_hier_n, &dd2->refine_z_hier));
1107: PetscCall(PetscArraycpy(dd2->refine_z_hier, dd->refine_z_hier, dd2->refine_z_hier_n));
1108: }
1109: if (dd->refine_y_hier) {
1110: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
1111: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
1112: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1113: PetscCall(PetscMalloc1(dd2->refine_y_hier_n, &dd2->refine_y_hier));
1114: PetscCall(PetscArraycpy(dd2->refine_y_hier, dd->refine_y_hier, dd2->refine_y_hier_n));
1115: }
1116: if (dd->refine_x_hier) {
1117: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
1118: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
1119: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1120: PetscCall(PetscMalloc1(dd2->refine_x_hier_n, &dd2->refine_x_hier));
1121: PetscCall(PetscArraycpy(dd2->refine_x_hier, dd->refine_x_hier, dd2->refine_x_hier_n));
1122: }
1124: /* copy vector type information */
1125: PetscCall(DMSetVecType(da2, da->vectype));
1127: dd2->lf = dd->lf;
1128: dd2->lj = dd->lj;
1130: da2->leveldown = da->leveldown;
1131: da2->levelup = da->levelup + 1;
1133: PetscCall(DMSetUp(da2));
1135: /* interpolate coordinates if they are set on the coarse grid */
1136: PetscCall(DMGetCoordinates(da, &coordsc));
1137: if (coordsc) {
1138: DM cdaf, cdac;
1139: Mat II;
1141: PetscCall(DMGetCoordinateDM(da, &cdac));
1142: PetscCall(DMGetCoordinateDM(da2, &cdaf));
1143: /* force creation of the coordinate vector */
1144: PetscCall(DMDASetUniformCoordinates(da2, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
1145: PetscCall(DMGetCoordinates(da2, &coordsf));
1146: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, NULL));
1147: PetscCall(MatInterpolate(II, coordsc, coordsf));
1148: PetscCall(MatDestroy(&II));
1149: }
1151: for (i = 0; i < da->bs; i++) {
1152: const char *fieldname;
1153: PetscCall(DMDAGetFieldName(da, i, &fieldname));
1154: PetscCall(DMDASetFieldName(da2, i, fieldname));
1155: }
1157: *daref = da2;
1158: PetscFunctionReturn(PETSC_SUCCESS);
1159: }
1161: PetscErrorCode DMCoarsen_DA(DM dmf, MPI_Comm comm, DM *dmc)
1162: {
1163: PetscInt M, N, P, i, dim;
1164: Vec coordsc, coordsf;
1165: DM dmc2;
1166: DM_DA *dd = (DM_DA *)dmf->data, *dd2;
1168: PetscFunctionBegin;
1170: PetscAssertPointer(dmc, 3);
1172: PetscCall(DMGetDimension(dmf, &dim));
1173: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1174: M = dd->M / dd->coarsen_x;
1175: } else {
1176: M = 1 + (dd->M - 1) / dd->coarsen_x;
1177: }
1178: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1179: if (dim > 1) {
1180: N = dd->N / dd->coarsen_y;
1181: } else {
1182: N = 1;
1183: }
1184: } else {
1185: N = 1 + (dd->N - 1) / dd->coarsen_y;
1186: }
1187: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1188: if (dim > 2) {
1189: P = dd->P / dd->coarsen_z;
1190: } else {
1191: P = 1;
1192: }
1193: } else {
1194: P = 1 + (dd->P - 1) / dd->coarsen_z;
1195: }
1196: PetscCall(DMDACreate(PetscObjectComm((PetscObject)dmf), &dmc2));
1197: PetscCall(DMSetOptionsPrefix(dmc2, ((PetscObject)dmf)->prefix));
1198: PetscCall(DMSetDimension(dmc2, dim));
1199: PetscCall(DMDASetSizes(dmc2, M, N, P));
1200: PetscCall(DMDASetNumProcs(dmc2, dd->m, dd->n, dd->p));
1201: PetscCall(DMDASetBoundaryType(dmc2, dd->bx, dd->by, dd->bz));
1202: PetscCall(DMDASetDof(dmc2, dd->w));
1203: PetscCall(DMDASetStencilType(dmc2, dd->stencil_type));
1204: PetscCall(DMDASetStencilWidth(dmc2, dd->s));
1205: if (dim == 3) {
1206: PetscInt *lx, *ly, *lz;
1207: PetscCall(PetscMalloc3(dd->m, &lx, dd->n, &ly, dd->p, &lz));
1208: PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx));
1209: PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_y, dd->n, dd->ly, ly));
1210: PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_z, dd->p, dd->lz, lz));
1211: PetscCall(DMDASetOwnershipRanges(dmc2, lx, ly, lz));
1212: PetscCall(PetscFree3(lx, ly, lz));
1213: } else if (dim == 2) {
1214: PetscInt *lx, *ly;
1215: PetscCall(PetscMalloc2(dd->m, &lx, dd->n, &ly));
1216: PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx));
1217: PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_y, dd->n, dd->ly, ly));
1218: PetscCall(DMDASetOwnershipRanges(dmc2, lx, ly, NULL));
1219: PetscCall(PetscFree2(lx, ly));
1220: } else if (dim == 1) {
1221: PetscInt *lx;
1222: PetscCall(PetscMalloc1(dd->m, &lx));
1223: PetscCall(DMDACoarsenOwnershipRanges(dmf, (PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0), dd->s, dd->coarsen_x, dd->m, dd->lx, lx));
1224: PetscCall(DMDASetOwnershipRanges(dmc2, lx, NULL, NULL));
1225: PetscCall(PetscFree(lx));
1226: }
1227: dd2 = (DM_DA *)dmc2->data;
1229: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1230: /* dmc2->ops->createinterpolation = dmf->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1231: dmc2->ops->creatematrix = dmf->ops->creatematrix;
1232: dmc2->ops->getcoloring = dmf->ops->getcoloring;
1233: dd2->interptype = dd->interptype;
1235: /* copy fill information if given */
1236: if (dd->dfill) {
1237: PetscCall(PetscMalloc1(dd->dfill[dd->w] + dd->w + 1, &dd2->dfill));
1238: PetscCall(PetscArraycpy(dd2->dfill, dd->dfill, dd->dfill[dd->w] + dd->w + 1));
1239: }
1240: if (dd->ofill) {
1241: PetscCall(PetscMalloc1(dd->ofill[dd->w] + dd->w + 1, &dd2->ofill));
1242: PetscCall(PetscArraycpy(dd2->ofill, dd->ofill, dd->ofill[dd->w] + dd->w + 1));
1243: }
1244: /* copy the refine information */
1245: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1246: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1247: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1249: if (dd->refine_z_hier) {
1250: if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_z_hier_n) dd2->refine_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 1];
1251: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_z_hier_n) dd2->coarsen_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 2];
1252: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1253: PetscCall(PetscMalloc1(dd2->refine_z_hier_n, &dd2->refine_z_hier));
1254: PetscCall(PetscArraycpy(dd2->refine_z_hier, dd->refine_z_hier, dd2->refine_z_hier_n));
1255: }
1256: if (dd->refine_y_hier) {
1257: if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_y_hier_n) dd2->refine_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 1];
1258: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_y_hier_n) dd2->coarsen_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 2];
1259: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1260: PetscCall(PetscMalloc1(dd2->refine_y_hier_n, &dd2->refine_y_hier));
1261: PetscCall(PetscArraycpy(dd2->refine_y_hier, dd->refine_y_hier, dd2->refine_y_hier_n));
1262: }
1263: if (dd->refine_x_hier) {
1264: if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_x_hier_n) dd2->refine_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 1];
1265: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_x_hier_n) dd2->coarsen_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 2];
1266: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1267: PetscCall(PetscMalloc1(dd2->refine_x_hier_n, &dd2->refine_x_hier));
1268: PetscCall(PetscArraycpy(dd2->refine_x_hier, dd->refine_x_hier, dd2->refine_x_hier_n));
1269: }
1271: /* copy vector type information */
1272: PetscCall(DMSetVecType(dmc2, dmf->vectype));
1274: dd2->lf = dd->lf;
1275: dd2->lj = dd->lj;
1277: dmc2->leveldown = dmf->leveldown + 1;
1278: dmc2->levelup = dmf->levelup;
1280: PetscCall(DMSetUp(dmc2));
1282: /* inject coordinates if they are set on the fine grid */
1283: PetscCall(DMGetCoordinates(dmf, &coordsf));
1284: if (coordsf) {
1285: DM cdaf, cdac;
1286: Mat inject;
1287: VecScatter vscat;
1289: PetscCall(DMGetCoordinateDM(dmf, &cdaf));
1290: PetscCall(DMGetCoordinateDM(dmc2, &cdac));
1291: /* force creation of the coordinate vector */
1292: PetscCall(DMDASetUniformCoordinates(dmc2, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
1293: PetscCall(DMGetCoordinates(dmc2, &coordsc));
1295: PetscCall(DMCreateInjection(cdac, cdaf, &inject));
1296: PetscCall(MatScatterGetVecScatter(inject, &vscat));
1297: PetscCall(VecScatterBegin(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD));
1298: PetscCall(VecScatterEnd(vscat, coordsf, coordsc, INSERT_VALUES, SCATTER_FORWARD));
1299: PetscCall(MatDestroy(&inject));
1300: }
1302: for (i = 0; i < dmf->bs; i++) {
1303: const char *fieldname;
1304: PetscCall(DMDAGetFieldName(dmf, i, &fieldname));
1305: PetscCall(DMDASetFieldName(dmc2, i, fieldname));
1306: }
1308: *dmc = dmc2;
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }
1312: PetscErrorCode DMRefineHierarchy_DA(DM da, PetscInt nlevels, DM daf[])
1313: {
1314: PetscInt i, n, *refx, *refy, *refz;
1316: PetscFunctionBegin;
1318: PetscCheck(nlevels >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "nlevels cannot be negative");
1319: if (nlevels == 0) PetscFunctionReturn(PETSC_SUCCESS);
1320: PetscAssertPointer(daf, 3);
1322: /* Get refinement factors, defaults taken from the coarse DMDA */
1323: PetscCall(PetscMalloc3(nlevels, &refx, nlevels, &refy, nlevels, &refz));
1324: for (i = 0; i < nlevels; i++) PetscCall(DMDAGetRefinementFactor(da, &refx[i], &refy[i], &refz[i]));
1325: n = nlevels;
1326: PetscCall(PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_x", refx, &n, NULL));
1327: n = nlevels;
1328: PetscCall(PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_y", refy, &n, NULL));
1329: n = nlevels;
1330: PetscCall(PetscOptionsGetIntArray(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_refine_hierarchy_z", refz, &n, NULL));
1332: PetscCall(DMDASetRefinementFactor(da, refx[0], refy[0], refz[0]));
1333: PetscCall(DMRefine(da, PetscObjectComm((PetscObject)da), &daf[0]));
1334: for (i = 1; i < nlevels; i++) {
1335: PetscCall(DMDASetRefinementFactor(daf[i - 1], refx[i], refy[i], refz[i]));
1336: PetscCall(DMRefine(daf[i - 1], PetscObjectComm((PetscObject)da), &daf[i]));
1337: }
1338: PetscCall(PetscFree3(refx, refy, refz));
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: PetscErrorCode DMCoarsenHierarchy_DA(DM da, PetscInt nlevels, DM dac[])
1343: {
1344: PetscInt i;
1346: PetscFunctionBegin;
1348: PetscCheck(nlevels >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "nlevels cannot be negative");
1349: if (nlevels == 0) PetscFunctionReturn(PETSC_SUCCESS);
1350: PetscAssertPointer(dac, 3);
1351: PetscCall(DMCoarsen(da, PetscObjectComm((PetscObject)da), &dac[0]));
1352: for (i = 1; i < nlevels; i++) PetscCall(DMCoarsen(dac[i - 1], PetscObjectComm((PetscObject)da), &dac[i]));
1353: PetscFunctionReturn(PETSC_SUCCESS);
1354: }
1356: static PetscErrorCode DMDASetGLLCoordinates_1d(DM dm, PetscInt n, PetscReal *nodes)
1357: {
1358: PetscInt i, j, xs, xn, q;
1359: PetscScalar *xx;
1360: PetscReal h;
1361: Vec x;
1362: DM_DA *da = (DM_DA *)dm->data;
1364: PetscFunctionBegin;
1365: if (da->bx != DM_BOUNDARY_PERIODIC) {
1366: PetscCall(DMDAGetInfo(dm, NULL, &q, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1367: q = (q - 1) / (n - 1); /* number of spectral elements */
1368: h = 2.0 / q;
1369: PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xn, NULL, NULL));
1370: xs = xs / (n - 1);
1371: xn = xn / (n - 1);
1372: PetscCall(DMDASetUniformCoordinates(dm, -1., 1., 0., 0., 0., 0.));
1373: PetscCall(DMGetCoordinates(dm, &x));
1374: PetscCall(DMDAVecGetArray(dm, x, &xx));
1376: /* loop over local spectral elements */
1377: for (j = xs; j < xs + xn; j++) {
1378: /*
1379: Except for the first process, each process starts on the second GLL point of the first element on that process
1380: */
1381: for (i = (j == xs && xs > 0) ? 1 : 0; i < n; i++) xx[j * (n - 1) + i] = -1.0 + h * j + h * (nodes[i] + 1.0) / 2.;
1382: }
1383: PetscCall(DMDAVecRestoreArray(dm, x, &xx));
1384: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not yet implemented for periodic");
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: /*@
1390: 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
1392: Collective
1394: Input Parameters:
1395: + da - the `DMDA` object
1396: . n - the number of GLL nodes
1397: - nodes - the GLL nodes
1399: Level: advanced
1401: Note:
1402: The parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1403: on each process much be divisible by the number of GLL elements needed per process. This depends on whether the `DMDA` is
1404: periodic or not.
1406: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDACreate()`, `PetscDTGaussLobattoLegendreQuadrature()`, `DMGetCoordinates()`
1407: @*/
1408: PetscErrorCode DMDASetGLLCoordinates(DM da, PetscInt n, PetscReal *nodes)
1409: {
1410: PetscFunctionBegin;
1411: if (da->dim == 1) {
1412: PetscCall(DMDASetGLLCoordinates_1d(da, n, nodes));
1413: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not yet implemented for 2 or 3d");
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: PetscErrorCode DMGetCompatibility_DA(DM da1, DM dm2, PetscBool *compatible, PetscBool *set)
1418: {
1419: DM_DA *dd1 = (DM_DA *)da1->data, *dd2;
1420: DM da2;
1421: DMType dmtype2;
1422: PetscBool isda, compatibleLocal;
1423: PetscInt i;
1425: PetscFunctionBegin;
1426: PetscCheck(da1->setupcalled, PetscObjectComm((PetscObject)da1), PETSC_ERR_ARG_WRONGSTATE, "DMSetUp() must be called on first DM before DMGetCompatibility()");
1427: PetscCall(DMGetType(dm2, &dmtype2));
1428: PetscCall(PetscStrcmp(dmtype2, DMDA, &isda));
1429: if (isda) {
1430: da2 = dm2;
1431: dd2 = (DM_DA *)da2->data;
1432: PetscCheck(da2->setupcalled, PetscObjectComm((PetscObject)da2), PETSC_ERR_ARG_WRONGSTATE, "DMSetUp() must be called on second DM before DMGetCompatibility()");
1433: compatibleLocal = (PetscBool)(da1->dim == da2->dim);
1434: if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
1435: /* Global size ranks Boundary type */
1436: if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
1437: if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
1438: if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
1439: if (compatibleLocal) {
1440: for (i = 0; i < dd1->m; ++i) { compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i])); /* Local size */ }
1441: }
1442: if (compatibleLocal && da1->dim > 1) {
1443: for (i = 0; i < dd1->n; ++i) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
1444: }
1445: if (compatibleLocal && da1->dim > 2) {
1446: for (i = 0; i < dd1->p; ++i) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
1447: }
1448: *compatible = compatibleLocal;
1449: *set = PETSC_TRUE;
1450: } else {
1451: /* Decline to determine compatibility with other DM types */
1452: *set = PETSC_FALSE;
1453: }
1454: PetscFunctionReturn(PETSC_SUCCESS);
1455: }