Actual source code: da.c
1: #include <petsc/private/dmdaimpl.h>
3: /*@
4: DMDASetSizes - Sets the number of grid points in the three dimensional directions
6: Logically Collective on da
8: Input Parameters:
9: + da - the DMDA
10: . M - the global X size
11: . N - the global Y size
12: - P - the global Z size
14: Level: intermediate
16: Developer Notes:
17: Since the dimension may not yet have been set the code cannot error check for non-positive Y and Z number of grid points
19: .seealso: PetscSplitOwnership()
20: @*/
21: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
22: {
23: DM_DA *dd = (DM_DA*)da->data;
34: dd->M = M;
35: dd->N = N;
36: dd->P = P;
37: return 0;
38: }
40: /*@
41: DMDASetNumProcs - Sets the number of processes in each dimension
43: Logically Collective on da
45: Input Parameters:
46: + da - the DMDA
47: . m - the number of X procs (or PETSC_DECIDE)
48: . n - the number of Y procs (or PETSC_DECIDE)
49: - p - the number of Z procs (or PETSC_DECIDE)
51: Level: intermediate
53: .seealso: DMDASetSizes(), PetscSplitOwnership()
54: @*/
55: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
56: {
57: DM_DA *dd = (DM_DA*)da->data;
64: dd->m = m;
65: dd->n = n;
66: dd->p = p;
67: if (da->dim == 2) {
68: PetscMPIInt size;
69: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
70: if ((dd->m > 0) && (dd->n < 0)) {
71: dd->n = size/dd->m;
73: }
74: if ((dd->n > 0) && (dd->m < 0)) {
75: dd->m = size/dd->n;
77: }
78: }
79: return 0;
80: }
82: /*@
83: DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.
85: Not collective
87: Input Parameters:
88: + da - The DMDA
89: - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC
91: Level: intermediate
93: .seealso: DMDACreate(), DMDestroy(), DMDA, DMBoundaryType
94: @*/
95: PetscErrorCode DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)
96: {
97: DM_DA *dd = (DM_DA*)da->data;
104: dd->bx = bx;
105: dd->by = by;
106: dd->bz = bz;
107: return 0;
108: }
110: /*@
111: DMDASetDof - Sets the number of degrees of freedom per vertex
113: Not collective
115: Input Parameters:
116: + da - The DMDA
117: - dof - Number of degrees of freedom
119: Level: intermediate
121: .seealso: DMDAGetDof(), DMDACreate(), DMDestroy(), DMDA
122: @*/
123: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
124: {
125: DM_DA *dd = (DM_DA*)da->data;
130: dd->w = dof;
131: da->bs = dof;
132: return 0;
133: }
135: /*@
136: DMDAGetDof - Gets the number of degrees of freedom per vertex
138: Not collective
140: Input Parameter:
141: . da - The DMDA
143: Output Parameter:
144: . dof - Number of degrees of freedom
146: Level: intermediate
148: .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
149: @*/
150: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
151: {
152: DM_DA *dd = (DM_DA *) da->data;
156: *dof = dd->w;
157: return 0;
158: }
160: /*@
161: DMDAGetOverlap - Gets the size of the per-processor overlap.
163: Not collective
165: Input Parameter:
166: . da - The DMDA
168: Output Parameters:
169: + x - Overlap in the x direction
170: . y - Overlap in the y direction
171: - z - Overlap in the z direction
173: Level: intermediate
175: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
176: @*/
177: PetscErrorCode DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
178: {
179: DM_DA *dd = (DM_DA*)da->data;
182: if (x) *x = dd->xol;
183: if (y) *y = dd->yol;
184: if (z) *z = dd->zol;
185: return 0;
186: }
188: /*@
189: DMDASetOverlap - Sets the size of the per-processor overlap.
191: Not collective
193: Input Parameters:
194: + da - The DMDA
195: . x - Overlap in the x direction
196: . y - Overlap in the y direction
197: - z - Overlap in the z direction
199: Level: intermediate
201: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
202: @*/
203: PetscErrorCode DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
204: {
205: DM_DA *dd = (DM_DA*)da->data;
211: dd->xol = x;
212: dd->yol = y;
213: dd->zol = z;
214: return 0;
215: }
217: /*@
218: DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.
220: Not collective
222: Input Parameters:
223: . da - The DMDA
225: Output Parameters:
226: . Nsub - Number of local subdomains created upon decomposition
228: Level: intermediate
230: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
231: @*/
232: PetscErrorCode DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
233: {
234: DM_DA *dd = (DM_DA*)da->data;
237: if (Nsub) *Nsub = dd->Nsub;
238: return 0;
239: }
241: /*@
242: DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.
244: Not collective
246: Input Parameters:
247: + da - The DMDA
248: - Nsub - The number of local subdomains requested
250: Level: intermediate
252: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
253: @*/
254: PetscErrorCode DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
255: {
256: DM_DA *dd = (DM_DA*)da->data;
260: dd->Nsub = Nsub;
261: return 0;
262: }
264: /*@
265: DMDASetOffset - Sets the index offset of the DA.
267: Collective on DA
269: Input Parameters:
270: + da - The DMDA
271: . xo - The offset in the x direction
272: . yo - The offset in the y direction
273: - zo - The offset in the z direction
275: Level: intermediate
277: Notes:
278: This is used primarily to overlap a computation on a local DA with that on a global DA without
279: changing boundary conditions or subdomain features that depend upon the global offsets.
281: .seealso: DMDAGetOffset(), DMDAVecGetArray()
282: @*/
283: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
284: {
285: DM_DA *dd = (DM_DA*)da->data;
294: dd->xo = xo;
295: dd->yo = yo;
296: dd->zo = zo;
297: dd->Mo = Mo;
298: dd->No = No;
299: dd->Po = Po;
301: if (da->coordinateDM) {
302: DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
303: }
304: return 0;
305: }
307: /*@
308: DMDAGetOffset - Gets the index offset of the DA.
310: Not collective
312: Input Parameter:
313: . da - The DMDA
315: Output Parameters:
316: + xo - The offset in the x direction
317: . yo - The offset in the y direction
318: . zo - The offset in the z direction
319: . Mo - The global size in the x direction
320: . No - The global size in the y direction
321: - Po - The global size in the z direction
323: Level: intermediate
325: .seealso: DMDASetOffset(), DMDAVecGetArray()
326: @*/
327: PetscErrorCode DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
328: {
329: DM_DA *dd = (DM_DA*)da->data;
332: if (xo) *xo = dd->xo;
333: if (yo) *yo = dd->yo;
334: if (zo) *zo = dd->zo;
335: if (Mo) *Mo = dd->Mo;
336: if (No) *No = dd->No;
337: if (Po) *Po = dd->Po;
338: return 0;
339: }
341: /*@
342: DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.
344: Not collective
346: Input Parameter:
347: . da - The DMDA
349: Output Parameters:
350: + xs - The start of the region in x
351: . ys - The start of the region in y
352: . zs - The start of the region in z
353: . xs - The size of the region in x
354: . ys - The size of the region in y
355: - zs - The size of the region in z
357: Level: intermediate
359: .seealso: DMDAGetOffset(), DMDAVecGetArray()
360: @*/
361: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
362: {
363: DM_DA *dd = (DM_DA*)da->data;
366: if (xs) *xs = dd->nonxs;
367: if (ys) *ys = dd->nonys;
368: if (zs) *zs = dd->nonzs;
369: if (xm) *xm = dd->nonxm;
370: if (ym) *ym = dd->nonym;
371: if (zm) *zm = dd->nonzm;
372: return 0;
373: }
375: /*@
376: DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.
378: Collective on DA
380: Input Parameters:
381: + da - The DMDA
382: . xs - The start of the region in x
383: . ys - The start of the region in y
384: . zs - The start of the region in z
385: . xs - The size of the region in x
386: . ys - The size of the region in y
387: - zs - The size of the region in z
389: Level: intermediate
391: .seealso: DMDAGetOffset(), DMDAVecGetArray()
392: @*/
393: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
394: {
395: DM_DA *dd = (DM_DA*)da->data;
404: dd->nonxs = xs;
405: dd->nonys = ys;
406: dd->nonzs = zs;
407: dd->nonxm = xm;
408: dd->nonym = ym;
409: dd->nonzm = zm;
411: return 0;
412: }
414: /*@
415: DMDASetStencilType - Sets the type of the communication stencil
417: Logically Collective on da
419: Input Parameters:
420: + da - The DMDA
421: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
423: Level: intermediate
425: .seealso: DMDACreate(), DMDestroy(), DMDA
426: @*/
427: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
428: {
429: DM_DA *dd = (DM_DA*)da->data;
434: dd->stencil_type = stype;
435: return 0;
436: }
438: /*@
439: DMDAGetStencilType - Gets the type of the communication stencil
441: Not collective
443: Input Parameter:
444: . da - The DMDA
446: Output Parameter:
447: . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
449: Level: intermediate
451: .seealso: DMDACreate(), DMDestroy(), DMDA
452: @*/
453: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
454: {
455: DM_DA *dd = (DM_DA*)da->data;
459: *stype = dd->stencil_type;
460: return 0;
461: }
463: /*@
464: DMDASetStencilWidth - Sets the width of the communication stencil
466: Logically Collective on da
468: Input Parameters:
469: + da - The DMDA
470: - width - The stencil width
472: Level: intermediate
474: .seealso: DMDACreate(), DMDestroy(), DMDA
475: @*/
476: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
477: {
478: DM_DA *dd = (DM_DA*)da->data;
483: dd->s = width;
484: return 0;
485: }
487: /*@
488: DMDAGetStencilWidth - Gets the width of the communication stencil
490: Not collective
492: Input Parameter:
493: . da - The DMDA
495: Output Parameter:
496: . width - The stencil width
498: Level: intermediate
500: .seealso: DMDACreate(), DMDestroy(), DMDA
501: @*/
502: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
503: {
504: DM_DA *dd = (DM_DA *) da->data;
508: *width = dd->s;
509: return 0;
510: }
512: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
513: {
514: PetscInt i,sum;
517: for (i=sum=0; i<m; i++) sum += lx[i];
519: return 0;
520: }
522: /*@
523: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
525: Logically Collective on da
527: Input Parameters:
528: + da - The DMDA
529: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
530: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
531: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.
533: Level: intermediate
535: Note: these numbers are NOT multiplied by the number of dof per node.
537: .seealso: DMDACreate(), DMDestroy(), DMDA
538: @*/
539: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
540: {
541: DM_DA *dd = (DM_DA*)da->data;
545: if (lx) {
547: DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
548: if (!dd->lx) {
549: PetscMalloc1(dd->m, &dd->lx);
550: }
551: PetscArraycpy(dd->lx, lx, dd->m);
552: }
553: if (ly) {
555: DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
556: if (!dd->ly) {
557: PetscMalloc1(dd->n, &dd->ly);
558: }
559: PetscArraycpy(dd->ly, ly, dd->n);
560: }
561: if (lz) {
563: DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
564: if (!dd->lz) {
565: PetscMalloc1(dd->p, &dd->lz);
566: }
567: PetscArraycpy(dd->lz, lz, dd->p);
568: }
569: return 0;
570: }
572: /*@
573: DMDASetInterpolationType - Sets the type of interpolation that will be
574: returned by DMCreateInterpolation()
576: Logically Collective on da
578: Input Parameters:
579: + da - initial distributed array
580: - ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
582: Level: intermediate
584: Notes:
585: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
587: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
588: @*/
589: PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
590: {
591: DM_DA *dd = (DM_DA*)da->data;
595: dd->interptype = ctype;
596: return 0;
597: }
599: /*@
600: DMDAGetInterpolationType - Gets the type of interpolation that will be
601: used by DMCreateInterpolation()
603: Not Collective
605: Input Parameter:
606: . da - distributed array
608: Output Parameter:
609: . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
611: Level: intermediate
613: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
614: @*/
615: PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
616: {
617: DM_DA *dd = (DM_DA*)da->data;
621: *ctype = dd->interptype;
622: return 0;
623: }
625: /*@C
626: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
627: processes neighbors.
629: Not Collective
631: Input Parameter:
632: . da - the DMDA object
634: Output Parameters:
635: . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
636: this process itself is in the list
638: Notes:
639: In 2d the array is of length 9, in 3d of length 27
640: Not supported in 1d
641: Do not free the array, it is freed when the DMDA is destroyed.
643: Fortran Notes:
644: In fortran you must pass in an array of the appropriate length.
646: Level: intermediate
648: @*/
649: PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
650: {
651: DM_DA *dd = (DM_DA*)da->data;
654: *ranks = dd->neighbors;
655: return 0;
656: }
658: /*@C
659: DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
661: Not Collective
663: Input Parameter:
664: . da - the DMDA object
666: Output Parameters:
667: + lx - ownership along x direction (optional)
668: . ly - ownership along y direction (optional)
669: - lz - ownership along z direction (optional)
671: Level: intermediate
673: Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
675: In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
676: eighth arguments from DMDAGetInfo()
678: In C you should not free these arrays, nor change the values in them. They will only have valid values while the
679: DMDA they came from still exists (has not been destroyed).
681: These numbers are NOT multiplied by the number of dof per node.
683: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
684: @*/
685: PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
686: {
687: DM_DA *dd = (DM_DA*)da->data;
690: if (lx) *lx = dd->lx;
691: if (ly) *ly = dd->ly;
692: if (lz) *lz = dd->lz;
693: return 0;
694: }
696: /*@
697: DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
699: Logically Collective on da
701: Input Parameters:
702: + da - the DMDA object
703: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
704: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
705: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
707: Options Database:
708: + -da_refine_x refine_x - refinement ratio in x direction
709: . -da_refine_y rafine_y - refinement ratio in y direction
710: . -da_refine_z refine_z - refinement ratio in z direction
711: - -da_refine <n> - refine the DMDA object n times when it is created.
713: Level: intermediate
715: Notes:
716: Pass PETSC_IGNORE to leave a value unchanged
718: .seealso: DMRefine(), DMDAGetRefinementFactor()
719: @*/
720: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
721: {
722: DM_DA *dd = (DM_DA*)da->data;
729: if (refine_x > 0) dd->refine_x = refine_x;
730: if (refine_y > 0) dd->refine_y = refine_y;
731: if (refine_z > 0) dd->refine_z = refine_z;
732: return 0;
733: }
735: /*@C
736: DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
738: Not Collective
740: Input Parameter:
741: . da - the DMDA object
743: Output Parameters:
744: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
745: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
746: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
748: Level: intermediate
750: Notes:
751: Pass NULL for values you do not need
753: .seealso: DMRefine(), DMDASetRefinementFactor()
754: @*/
755: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
756: {
757: DM_DA *dd = (DM_DA*)da->data;
760: if (refine_x) *refine_x = dd->refine_x;
761: if (refine_y) *refine_y = dd->refine_y;
762: if (refine_z) *refine_z = dd->refine_z;
763: return 0;
764: }
766: /*@C
767: DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
769: Logically Collective on da
771: Input Parameters:
772: + da - the DMDA object
773: - f - the function that allocates the matrix for that specific DMDA
775: Level: developer
777: Notes:
778: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
779: the diagonal and off-diagonal blocks of the matrix
781: Not supported from Fortran
783: .seealso: DMCreateMatrix(), DMDASetBlockFills()
784: @*/
785: PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
786: {
788: da->ops->creatematrix = f;
789: return 0;
790: }
792: /*
793: Creates "balanced" ownership ranges after refinement, constrained by the need for the
794: fine grid boundaries to fall within one stencil width of the coarse partition.
796: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
797: */
798: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
799: {
800: PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
803: if (ratio == 1) {
804: PetscArraycpy(lf,lc,m);
805: return 0;
806: }
807: for (i=0; i<m; i++) totalc += lc[i];
808: remaining = (!periodic) + ratio * (totalc - (!periodic));
809: for (i=0; i<m; i++) {
810: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
811: if (i == m-1) lf[i] = want;
812: else {
813: const PetscInt nextc = startc + lc[i];
814: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
815: * coarse stencil width of the first coarse node in the next subdomain. */
816: while ((startf+want)/ratio < nextc - stencil_width) want++;
817: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
818: * coarse stencil width of the last coarse node in the current subdomain. */
819: while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
820: /* Make sure all constraints are satisfied */
821: if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
822: || ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible refined ownership range");
823: }
824: lf[i] = want;
825: startc += lc[i];
826: startf += lf[i];
827: remaining -= lf[i];
828: }
829: return 0;
830: }
832: /*
833: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
834: fine grid boundaries to fall within one stencil width of the coarse partition.
836: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
837: */
838: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
839: {
840: PetscInt i,totalf,remaining,startc,startf;
843: if (ratio == 1) {
844: PetscArraycpy(lc,lf,m);
845: return 0;
846: }
847: for (i=0,totalf=0; i<m; i++) totalf += lf[i];
848: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
849: for (i=0,startc=0,startf=0; i<m; i++) {
850: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
851: if (i == m-1) lc[i] = want;
852: else {
853: const PetscInt nextf = startf+lf[i];
854: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
855: * node is within one stencil width. */
856: while (nextf/ratio < startc+want-stencil_width) want--;
857: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
858: * fine node is within one stencil width. */
859: while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
860: if (want < 0 || want > remaining
861: || (nextf/ratio < startc+want-stencil_width) || ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible coarsened ownership range");
862: }
863: lc[i] = want;
864: startc += lc[i];
865: startf += lf[i];
866: remaining -= lc[i];
867: }
868: return 0;
869: }
871: PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
872: {
873: PetscInt M,N,P,i,dim;
874: DM da2;
875: DM_DA *dd = (DM_DA*)da->data,*dd2;
880: DMGetDimension(da, &dim);
881: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
882: M = dd->refine_x*dd->M;
883: } else {
884: M = 1 + dd->refine_x*(dd->M - 1);
885: }
886: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
887: if (dim > 1) {
888: N = dd->refine_y*dd->N;
889: } else {
890: N = 1;
891: }
892: } else {
893: N = 1 + dd->refine_y*(dd->N - 1);
894: }
895: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
896: if (dim > 2) {
897: P = dd->refine_z*dd->P;
898: } else {
899: P = 1;
900: }
901: } else {
902: P = 1 + dd->refine_z*(dd->P - 1);
903: }
904: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
905: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
906: DMSetDimension(da2,dim);
907: DMDASetSizes(da2,M,N,P);
908: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
909: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
910: DMDASetDof(da2,dd->w);
911: DMDASetStencilType(da2,dd->stencil_type);
912: DMDASetStencilWidth(da2,dd->s);
913: if (dim == 3) {
914: PetscInt *lx,*ly,*lz;
915: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
916: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
917: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
918: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
919: DMDASetOwnershipRanges(da2,lx,ly,lz);
920: PetscFree3(lx,ly,lz);
921: } else if (dim == 2) {
922: PetscInt *lx,*ly;
923: PetscMalloc2(dd->m,&lx,dd->n,&ly);
924: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
925: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
926: DMDASetOwnershipRanges(da2,lx,ly,NULL);
927: PetscFree2(lx,ly);
928: } else if (dim == 1) {
929: PetscInt *lx;
930: PetscMalloc1(dd->m,&lx);
931: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
932: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
933: PetscFree(lx);
934: }
935: dd2 = (DM_DA*)da2->data;
937: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
938: da2->ops->creatematrix = da->ops->creatematrix;
939: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
940: da2->ops->getcoloring = da->ops->getcoloring;
941: dd2->interptype = dd->interptype;
943: /* copy fill information if given */
944: if (dd->dfill) {
945: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
946: PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);
947: }
948: if (dd->ofill) {
949: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
950: PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);
951: }
952: /* copy the refine information */
953: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
954: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
955: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
957: if (dd->refine_z_hier) {
958: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) {
959: dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
960: }
961: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) {
962: dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
963: }
964: dd2->refine_z_hier_n = dd->refine_z_hier_n;
965: PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
966: PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);
967: }
968: if (dd->refine_y_hier) {
969: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) {
970: dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
971: }
972: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) {
973: dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
974: }
975: dd2->refine_y_hier_n = dd->refine_y_hier_n;
976: PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
977: PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);
978: }
979: if (dd->refine_x_hier) {
980: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) {
981: dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
982: }
983: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) {
984: dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
985: }
986: dd2->refine_x_hier_n = dd->refine_x_hier_n;
987: PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
988: PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);
989: }
991: /* copy vector type information */
992: DMSetVecType(da2,da->vectype);
994: dd2->lf = dd->lf;
995: dd2->lj = dd->lj;
997: da2->leveldown = da->leveldown;
998: da2->levelup = da->levelup + 1;
1000: DMSetUp(da2);
1002: /* interpolate coordinates if they are set on the coarse grid */
1003: if (da->coordinates) {
1004: DM cdaf,cdac;
1005: Vec coordsc,coordsf;
1006: Mat II;
1008: DMGetCoordinateDM(da,&cdac);
1009: DMGetCoordinates(da,&coordsc);
1010: DMGetCoordinateDM(da2,&cdaf);
1011: /* force creation of the coordinate vector */
1012: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1013: DMGetCoordinates(da2,&coordsf);
1014: DMCreateInterpolation(cdac,cdaf,&II,NULL);
1015: MatInterpolate(II,coordsc,coordsf);
1016: MatDestroy(&II);
1017: }
1019: for (i=0; i<da->bs; i++) {
1020: const char *fieldname;
1021: DMDAGetFieldName(da,i,&fieldname);
1022: DMDASetFieldName(da2,i,fieldname);
1023: }
1025: *daref = da2;
1026: return 0;
1027: }
1029: PetscErrorCode DMCoarsen_DA(DM dmf, MPI_Comm comm,DM *dmc)
1030: {
1031: PetscInt M,N,P,i,dim;
1032: DM dmc2;
1033: DM_DA *dd = (DM_DA*)dmf->data,*dd2;
1038: DMGetDimension(dmf, &dim);
1039: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1040: M = dd->M / dd->coarsen_x;
1041: } else {
1042: M = 1 + (dd->M - 1) / dd->coarsen_x;
1043: }
1044: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1045: if (dim > 1) {
1046: N = dd->N / dd->coarsen_y;
1047: } else {
1048: N = 1;
1049: }
1050: } else {
1051: N = 1 + (dd->N - 1) / dd->coarsen_y;
1052: }
1053: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1054: if (dim > 2) {
1055: P = dd->P / dd->coarsen_z;
1056: } else {
1057: P = 1;
1058: }
1059: } else {
1060: P = 1 + (dd->P - 1) / dd->coarsen_z;
1061: }
1062: DMDACreate(PetscObjectComm((PetscObject)dmf),&dmc2);
1063: DMSetOptionsPrefix(dmc2,((PetscObject)dmf)->prefix);
1064: DMSetDimension(dmc2,dim);
1065: DMDASetSizes(dmc2,M,N,P);
1066: DMDASetNumProcs(dmc2,dd->m,dd->n,dd->p);
1067: DMDASetBoundaryType(dmc2,dd->bx,dd->by,dd->bz);
1068: DMDASetDof(dmc2,dd->w);
1069: DMDASetStencilType(dmc2,dd->stencil_type);
1070: DMDASetStencilWidth(dmc2,dd->s);
1071: if (dim == 3) {
1072: PetscInt *lx,*ly,*lz;
1073: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1074: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1075: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1076: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1077: DMDASetOwnershipRanges(dmc2,lx,ly,lz);
1078: PetscFree3(lx,ly,lz);
1079: } else if (dim == 2) {
1080: PetscInt *lx,*ly;
1081: PetscMalloc2(dd->m,&lx,dd->n,&ly);
1082: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1083: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1084: DMDASetOwnershipRanges(dmc2,lx,ly,NULL);
1085: PetscFree2(lx,ly);
1086: } else if (dim == 1) {
1087: PetscInt *lx;
1088: PetscMalloc1(dd->m,&lx);
1089: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1090: DMDASetOwnershipRanges(dmc2,lx,NULL,NULL);
1091: PetscFree(lx);
1092: }
1093: dd2 = (DM_DA*)dmc2->data;
1095: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1096: /* dmc2->ops->createinterpolation = dmf->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1097: dmc2->ops->creatematrix = dmf->ops->creatematrix;
1098: dmc2->ops->getcoloring = dmf->ops->getcoloring;
1099: dd2->interptype = dd->interptype;
1101: /* copy fill information if given */
1102: if (dd->dfill) {
1103: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1104: PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);
1105: }
1106: if (dd->ofill) {
1107: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1108: PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);
1109: }
1110: /* copy the refine information */
1111: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1112: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1113: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1115: if (dd->refine_z_hier) {
1116: if (dmf->levelup - dmf->leveldown -1 > -1 && dmf->levelup - dmf->leveldown - 1< dd->refine_z_hier_n) {
1117: dd2->refine_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 1];
1118: }
1119: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_z_hier_n) {
1120: dd2->coarsen_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 2];
1121: }
1122: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1123: PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1124: PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);
1125: }
1126: if (dd->refine_y_hier) {
1127: if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1< dd->refine_y_hier_n) {
1128: dd2->refine_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 1];
1129: }
1130: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_y_hier_n) {
1131: dd2->coarsen_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 2];
1132: }
1133: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1134: PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1135: PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);
1136: }
1137: if (dd->refine_x_hier) {
1138: if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_x_hier_n) {
1139: dd2->refine_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 1];
1140: }
1141: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_x_hier_n) {
1142: dd2->coarsen_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 2];
1143: }
1144: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1145: PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1146: PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);
1147: }
1149: /* copy vector type information */
1150: DMSetVecType(dmc2,dmf->vectype);
1152: dd2->lf = dd->lf;
1153: dd2->lj = dd->lj;
1155: dmc2->leveldown = dmf->leveldown + 1;
1156: dmc2->levelup = dmf->levelup;
1158: DMSetUp(dmc2);
1160: /* inject coordinates if they are set on the fine grid */
1161: if (dmf->coordinates) {
1162: DM cdaf,cdac;
1163: Vec coordsc,coordsf;
1164: Mat inject;
1165: VecScatter vscat;
1167: DMGetCoordinateDM(dmf,&cdaf);
1168: DMGetCoordinates(dmf,&coordsf);
1169: DMGetCoordinateDM(dmc2,&cdac);
1170: /* force creation of the coordinate vector */
1171: DMDASetUniformCoordinates(dmc2,0.0,1.0,0.0,1.0,0.0,1.0);
1172: DMGetCoordinates(dmc2,&coordsc);
1174: DMCreateInjection(cdac,cdaf,&inject);
1175: MatScatterGetVecScatter(inject,&vscat);
1176: VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1177: VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1178: MatDestroy(&inject);
1179: }
1181: for (i=0; i<dmf->bs; i++) {
1182: const char *fieldname;
1183: DMDAGetFieldName(dmf,i,&fieldname);
1184: DMDASetFieldName(dmc2,i,fieldname);
1185: }
1187: *dmc = dmc2;
1188: return 0;
1189: }
1191: PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1192: {
1193: PetscInt i,n,*refx,*refy,*refz;
1197: if (nlevels == 0) return 0;
1200: /* Get refinement factors, defaults taken from the coarse DMDA */
1201: PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1202: for (i=0; i<nlevels; i++) {
1203: DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1204: }
1205: n = nlevels;
1206: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1207: n = nlevels;
1208: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1209: n = nlevels;
1210: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);
1212: DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1213: DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1214: for (i=1; i<nlevels; i++) {
1215: DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1216: DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1217: }
1218: PetscFree3(refx,refy,refz);
1219: return 0;
1220: }
1222: PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1223: {
1224: PetscInt i;
1228: if (nlevels == 0) return 0;
1230: DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1231: for (i=1; i<nlevels; i++) {
1232: DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1233: }
1234: return 0;
1235: }
1237: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm,PetscInt n,PetscReal *nodes)
1238: {
1239: PetscInt i,j,xs,xn,q;
1240: PetscScalar *xx;
1241: PetscReal h;
1242: Vec x;
1243: DM_DA *da = (DM_DA*)dm->data;
1245: if (da->bx != DM_BOUNDARY_PERIODIC) {
1246: DMDAGetInfo(dm,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1247: q = (q-1)/(n-1); /* number of spectral elements */
1248: h = 2.0/q;
1249: DMDAGetCorners(dm,&xs,NULL,NULL,&xn,NULL,NULL);
1250: xs = xs/(n-1);
1251: xn = xn/(n-1);
1252: DMDASetUniformCoordinates(dm,-1.,1.,0.,0.,0.,0.);
1253: DMGetCoordinates(dm,&x);
1254: DMDAVecGetArray(dm,x,&xx);
1256: /* loop over local spectral elements */
1257: for (j=xs; j<xs+xn; j++) {
1258: /*
1259: Except for the first process, each process starts on the second GLL point of the first element on that process
1260: */
1261: for (i= (j == xs && xs > 0)? 1 : 0; i<n; i++) {
1262: xx[j*(n-1) + i] = -1.0 + h*j + h*(nodes[i]+1.0)/2.;
1263: }
1264: }
1265: DMDAVecRestoreArray(dm,x,&xx);
1266: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for periodic");
1267: return 0;
1268: }
1270: /*@
1272: 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
1274: Collective on da
1276: Input Parameters:
1277: + da - the DMDA object
1278: - n - the number of GLL nodes
1279: - nodes - the GLL nodes
1281: Notes:
1282: the parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1283: on each process much be divisible by the number of GLL elements needed per process. This depends on whether the DM is
1284: periodic or not.
1286: Level: advanced
1288: .seealso: DMDACreate(), PetscDTGaussLobattoLegendreQuadrature(), DMGetCoordinates()
1289: @*/
1290: PetscErrorCode DMDASetGLLCoordinates(DM da,PetscInt n,PetscReal *nodes)
1291: {
1292: if (da->dim == 1) {
1293: DMDASetGLLCoordinates_1d(da,n,nodes);
1294: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for 2 or 3d");
1295: return 0;
1296: }
1298: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM da1,DM dm2,PetscBool *compatible,PetscBool *set)
1299: {
1300: DM_DA *dd1 = (DM_DA*)da1->data,*dd2;
1301: DM da2;
1302: DMType dmtype2;
1303: PetscBool isda,compatibleLocal;
1304: PetscInt i;
1307: DMGetType(dm2,&dmtype2);
1308: PetscStrcmp(dmtype2,DMDA,&isda);
1309: if (isda) {
1310: da2 = dm2;
1311: dd2 = (DM_DA*)da2->data;
1313: compatibleLocal = (PetscBool)(da1->dim == da2->dim);
1314: if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
1315: /* Global size ranks Boundary type */
1316: if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
1317: if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
1318: if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
1319: if (compatibleLocal) {
1320: for (i=0; i<dd1->m; ++i) {
1321: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i])); /* Local size */
1322: }
1323: }
1324: if (compatibleLocal && da1->dim > 1) {
1325: for (i=0; i<dd1->n; ++i) {
1326: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
1327: }
1328: }
1329: if (compatibleLocal && da1->dim > 2) {
1330: for (i=0; i<dd1->p; ++i) {
1331: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
1332: }
1333: }
1334: *compatible = compatibleLocal;
1335: *set = PETSC_TRUE;
1336: } else {
1337: /* Decline to determine compatibility with other DM types */
1338: *set = PETSC_FALSE;
1339: }
1340: return 0;
1341: }