Actual source code: da.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/
6: /*@
7: DMDASetDim - Sets the dimension
9: Collective on DMDA
11: Input Parameters:
12: + da - the DMDA
13: - dim - the dimension (or PETSC_DECIDE)
15: Level: intermediate
17: .seealso: DaGetDim(), DMDASetSizes()
18: @*/
19: PetscErrorCode DMDASetDim(DM da, PetscInt dim)
20: {
21: DM_DA *dd = (DM_DA*)da->data;
25: if (dd->dim > 0 && dim != dd->dim) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot change DMDA dim from %D after it was set to %D",dd->dim,dim);
26: dd->dim = dim;
27: return(0);
28: }
32: /*@
33: DMDASetSizes - Sets the global sizes
35: Logically Collective on DMDA
37: Input Parameters:
38: + da - the DMDA
39: . M - the global X size (or PETSC_DECIDE)
40: . N - the global Y size (or PETSC_DECIDE)
41: - P - the global Z size (or PETSC_DECIDE)
43: Level: intermediate
45: .seealso: DMDAGetSize(), PetscSplitOwnership()
46: @*/
47: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
48: {
49: DM_DA *dd = (DM_DA*)da->data;
56: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
58: dd->M = M;
59: dd->N = N;
60: dd->P = P;
61: return(0);
62: }
66: /*@
67: DMDASetNumProcs - Sets the number of processes in each dimension
69: Logically Collective on DMDA
71: Input Parameters:
72: + da - the DMDA
73: . m - the number of X procs (or PETSC_DECIDE)
74: . n - the number of Y procs (or PETSC_DECIDE)
75: - p - the number of Z procs (or PETSC_DECIDE)
77: Level: intermediate
79: .seealso: DMDASetSizes(), DMDAGetSize(), PetscSplitOwnership()
80: @*/
81: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
82: {
83: DM_DA *dd = (DM_DA*)da->data;
91: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
92: dd->m = m;
93: dd->n = n;
94: dd->p = p;
95: if (dd->dim == 2) {
96: PetscMPIInt size;
97: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
98: if ((dd->m > 0) && (dd->n < 0)) {
99: dd->n = size/dd->m;
100: if (dd->n*dd->m != size) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"%D processes in X direction not divisible into comm size %d",m,size);
101: }
102: if ((dd->n > 0) && (dd->m < 0)) {
103: dd->m = size/dd->n;
104: if (dd->n*dd->m != size) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"%D processes in Y direction not divisible into comm size %d",n,size);
105: }
106: }
107: return(0);
108: }
112: /*@
113: DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.
115: Not collective
117: Input Parameter:
118: + da - The DMDA
119: - bx,by,bz - One of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC
121: Level: intermediate
123: .keywords: distributed array, periodicity
124: .seealso: DMDACreate(), DMDestroy(), DMDA, DMDABoundaryType
125: @*/
126: PetscErrorCode DMDASetBoundaryType(DM da,DMDABoundaryType bx,DMDABoundaryType by,DMDABoundaryType bz)
127: {
128: DM_DA *dd = (DM_DA*)da->data;
135: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
136: dd->bx = bx;
137: dd->by = by;
138: dd->bz = bz;
139: return(0);
140: }
144: /*@
145: DMDASetDof - Sets the number of degrees of freedom per vertex
147: Not collective
149: Input Parameter:
150: + da - The DMDA
151: - dof - Number of degrees of freedom
153: Level: intermediate
155: .keywords: distributed array, degrees of freedom
156: .seealso: DMDACreate(), DMDestroy(), DMDA
157: @*/
158: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
159: {
160: DM_DA *dd = (DM_DA*)da->data;
165: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
166: dd->w = dof;
167: da->bs = dof;
168: return(0);
169: }
173: /*@
174: DMDAGetOverlap - Gets the size of the per-processor overlap.
176: Not collective
178: Input Parameters:
179: . da - The DMDA
181: Output Parameters:
182: + x - Overlap in the x direction
183: . y - Overlap in the y direction
184: - z - Overlap in the z direction
186: Level: intermediate
188: .keywords: distributed array, overlap, domain decomposition
189: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
190: @*/
191: PetscErrorCode DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
192: {
193: DM_DA *dd = (DM_DA*)da->data;
197: if (x) *x = dd->xol;
198: if (y) *y = dd->yol;
199: if (z) *z = dd->zol;
200: return(0);
201: }
205: /*@
206: DMDASetOverlap - Sets the size of the per-processor overlap.
208: Not collective
210: Input Parameters:
211: + da - The DMDA
212: . x - Overlap in the x direction
213: . y - Overlap in the y direction
214: - z - Overlap in the z direction
216: Level: intermediate
218: .keywords: distributed array, overlap, domain decomposition
219: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
220: @*/
221: PetscErrorCode DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
222: {
223: DM_DA *dd = (DM_DA*)da->data;
230: dd->xol = x;
231: dd->yol = y;
232: dd->zol = z;
233: return(0);
234: }
239: /*@
240: DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.
242: Not collective
244: Input Parameters:
245: . da - The DMDA
247: Output Parameters:
248: + Nsub - Number of local subdomains created upon decomposition
250: Level: intermediate
252: .keywords: distributed array, domain decomposition
253: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
254: @*/
255: PetscErrorCode DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
256: {
257: DM_DA *dd = (DM_DA*)da->data;
261: if (Nsub) *Nsub = dd->Nsub;
262: return(0);
263: }
267: /*@
268: DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.
270: Not collective
272: Input Parameters:
273: + da - The DMDA
274: - Nsub - The number of local subdomains requested
276: Level: intermediate
278: .keywords: distributed array, domain decomposition
279: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
280: @*/
281: PetscErrorCode DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
282: {
283: DM_DA *dd = (DM_DA*)da->data;
288: dd->Nsub = Nsub;
289: return(0);
290: }
294: /*@
295: DMDASetOffset - Sets the index offset of the DA.
297: Collective on DA
299: Input Parameter:
300: + da - The DMDA
301: . xo - The offset in the x direction
302: . yo - The offset in the y direction
303: - zo - The offset in the z direction
305: Level: intermediate
307: Notes: This is used primarily to overlap a computation on a local DA with that on a global DA without
308: changing boundary conditions or subdomain features that depend upon the global offsets.
310: .keywords: distributed array, degrees of freedom
311: .seealso: DMDAGetOffset(), DMDAVecGetArray()
312: @*/
313: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
314: {
316: DM_DA *dd = (DM_DA*)da->data;
326: dd->xo = xo;
327: dd->yo = yo;
328: dd->zo = zo;
329: dd->Mo = Mo;
330: dd->No = No;
331: dd->Po = Po;
333: if (da->coordinateDM) {
334: DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
335: }
336: return(0);
337: }
341: /*@
342: DMDAGetOffset - Gets the index offset of the DA.
344: Not collective
346: Input Parameter:
347: . da - The DMDA
349: Output Parameters:
350: + xo - The offset in the x direction
351: . yo - The offset in the y direction
352: . zo - The offset in the z direction
353: . Mo - The global size in the x direction
354: . No - The global size in the y direction
355: - Po - The global size in the z direction
357: Level: intermediate
359: .keywords: distributed array, degrees of freedom
360: .seealso: DMDASetOffset(), DMDAVecGetArray()
361: @*/
362: PetscErrorCode DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
363: {
364: DM_DA *dd = (DM_DA*)da->data;
368: if (xo) *xo = dd->xo;
369: if (yo) *yo = dd->yo;
370: if (zo) *zo = dd->zo;
371: if (Mo) *Mo = dd->Mo;
372: if (No) *No = dd->No;
373: if (Po) *Po = dd->Po;
374: return(0);
375: }
379: /*@
380: DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.
382: Not collective
384: Input Parameter:
385: . da - The DMDA
387: Output Parameters:
388: + xs - The start of the region in x
389: . ys - The start of the region in y
390: . zs - The start of the region in z
391: . xs - The size of the region in x
392: . ys - The size of the region in y
393: . zs - The size of the region in z
395: Level: intermediate
397: .keywords: distributed array, degrees of freedom
398: .seealso: DMDAGetOffset(), DMDAVecGetArray()
399: @*/
400: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
401: {
402: DM_DA *dd = (DM_DA*)da->data;
406: if (xs) *xs = dd->nonxs;
407: if (ys) *ys = dd->nonys;
408: if (zs) *zs = dd->nonzs;
409: if (xm) *xm = dd->nonxm;
410: if (ym) *ym = dd->nonym;
411: if (zm) *zm = dd->nonzm;
412: return(0);
413: }
418: /*@
419: DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.
421: Collective on DA
423: Input Parameter:
424: + da - The DMDA
425: . xs - The start of the region in x
426: . ys - The start of the region in y
427: . zs - The start of the region in z
428: . xs - The size of the region in x
429: . ys - The size of the region in y
430: . zs - The size of the region in z
432: Level: intermediate
434: .keywords: distributed array, degrees of freedom
435: .seealso: DMDAGetOffset(), DMDAVecGetArray()
436: @*/
437: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
438: {
439: DM_DA *dd = (DM_DA*)da->data;
449: dd->nonxs = xs;
450: dd->nonys = ys;
451: dd->nonzs = zs;
452: dd->nonxm = xm;
453: dd->nonym = ym;
454: dd->nonzm = zm;
456: return(0);
457: }
461: /*@
462: DMDASetStencilType - Sets the type of the communication stencil
464: Logically Collective on DMDA
466: Input Parameter:
467: + da - The DMDA
468: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
470: Level: intermediate
472: .keywords: distributed array, stencil
473: .seealso: DMDACreate(), DMDestroy(), DMDA
474: @*/
475: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
476: {
477: DM_DA *dd = (DM_DA*)da->data;
482: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
483: dd->stencil_type = stype;
484: return(0);
485: }
489: /*@
490: DMDASetStencilWidth - Sets the width of the communication stencil
492: Logically Collective on DMDA
494: Input Parameter:
495: + da - The DMDA
496: - width - The stencil width
498: Level: intermediate
500: .keywords: distributed array, stencil
501: .seealso: DMDACreate(), DMDestroy(), DMDA
502: @*/
503: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
504: {
505: DM_DA *dd = (DM_DA*)da->data;
510: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
511: dd->s = width;
512: return(0);
513: }
517: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
518: {
519: PetscInt i,sum;
522: if (M < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
523: for (i=sum=0; i<m; i++) sum += lx[i];
524: if (sum != M) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
525: return(0);
526: }
530: /*@
531: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
533: Logically Collective on DMDA
535: Input Parameter:
536: + da - The DMDA
537: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
538: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
539: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.
541: Level: intermediate
543: Note: these numbers are NOT multiplied by the number of dof per node.
545: .keywords: distributed array
546: .seealso: DMDACreate(), DMDestroy(), DMDA
547: @*/
548: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
549: {
551: DM_DA *dd = (DM_DA*)da->data;
555: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
556: if (lx) {
557: if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
558: DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
559: if (!dd->lx) {
560: PetscMalloc(dd->m*sizeof(PetscInt), &dd->lx);
561: }
562: PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
563: }
564: if (ly) {
565: if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
566: DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
567: if (!dd->ly) {
568: PetscMalloc(dd->n*sizeof(PetscInt), &dd->ly);
569: }
570: PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
571: }
572: if (lz) {
573: if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
574: DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
575: if (!dd->lz) {
576: PetscMalloc(dd->p*sizeof(PetscInt), &dd->lz);
577: }
578: PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
579: }
580: return(0);
581: }
585: /*@
586: DMDASetInterpolationType - Sets the type of interpolation that will be
587: returned by DMCreateInterpolation()
589: Logically Collective on DMDA
591: Input Parameter:
592: + da - initial distributed array
593: . ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
595: Level: intermediate
597: Notes: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
599: .keywords: distributed array, interpolation
601: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
602: @*/
603: PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
604: {
605: DM_DA *dd = (DM_DA*)da->data;
610: dd->interptype = ctype;
611: return(0);
612: }
616: /*@
617: DMDAGetInterpolationType - Gets the type of interpolation that will be
618: used by DMCreateInterpolation()
620: Not Collective
622: Input Parameter:
623: . da - distributed array
625: Output Parameter:
626: . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
628: Level: intermediate
630: .keywords: distributed array, interpolation
632: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
633: @*/
634: PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
635: {
636: DM_DA *dd = (DM_DA*)da->data;
641: *ctype = dd->interptype;
642: return(0);
643: }
647: /*@C
648: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
649: processes neighbors.
651: Not Collective
653: Input Parameter:
654: . da - the DMDA object
656: Output Parameters:
657: . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
658: this process itself is in the list
660: Notes: In 2d the array is of length 9, in 3d of length 27
661: Not supported in 1d
662: Do not free the array, it is freed when the DMDA is destroyed.
664: Fortran Notes: In fortran you must pass in an array of the appropriate length.
666: Level: intermediate
668: @*/
669: PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
670: {
671: DM_DA *dd = (DM_DA*)da->data;
675: *ranks = dd->neighbors;
676: return(0);
677: }
681: /*@C
682: DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
684: Not Collective
686: Input Parameter:
687: . da - the DMDA object
689: Output Parameter:
690: + lx - ownership along x direction (optional)
691: . ly - ownership along y direction (optional)
692: - lz - ownership along z direction (optional)
694: Level: intermediate
696: Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
698: In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
699: eighth arguments from DMDAGetInfo()
701: In C you should not free these arrays, nor change the values in them. They will only have valid values while the
702: DMDA they came from still exists (has not been destroyed).
704: These numbers are NOT multiplied by the number of dof per node.
706: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
707: @*/
708: PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
709: {
710: DM_DA *dd = (DM_DA*)da->data;
714: if (lx) *lx = dd->lx;
715: if (ly) *ly = dd->ly;
716: if (lz) *lz = dd->lz;
717: return(0);
718: }
722: /*@
723: DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
725: Logically Collective on DMDA
727: Input Parameters:
728: + da - the DMDA object
729: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
730: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
731: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
733: Options Database:
734: + -da_refine_x - refinement ratio in x direction
735: . -da_refine_y - refinement ratio in y direction
736: - -da_refine_z - refinement ratio in z direction
738: Level: intermediate
740: Notes: Pass PETSC_IGNORE to leave a value unchanged
742: .seealso: DMRefine(), DMDAGetRefinementFactor()
743: @*/
744: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
745: {
746: DM_DA *dd = (DM_DA*)da->data;
754: if (refine_x > 0) dd->refine_x = refine_x;
755: if (refine_y > 0) dd->refine_y = refine_y;
756: if (refine_z > 0) dd->refine_z = refine_z;
757: return(0);
758: }
762: /*@C
763: DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
765: Not Collective
767: Input Parameter:
768: . da - the DMDA object
770: Output Parameters:
771: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
772: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
773: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
775: Level: intermediate
777: Notes: Pass NULL for values you do not need
779: .seealso: DMRefine(), DMDASetRefinementFactor()
780: @*/
781: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
782: {
783: DM_DA *dd = (DM_DA*)da->data;
787: if (refine_x) *refine_x = dd->refine_x;
788: if (refine_y) *refine_y = dd->refine_y;
789: if (refine_z) *refine_z = dd->refine_z;
790: return(0);
791: }
795: /*@C
796: DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
798: Logically Collective on DMDA
800: Input Parameters:
801: + da - the DMDA object
802: - f - the function that allocates the matrix for that specific DMDA
804: Level: developer
806: Notes: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
807: the diagonal and off-diagonal blocks of the matrix
809: .seealso: DMCreateMatrix(), DMDASetBlockFills()
810: @*/
811: PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, MatType,Mat*))
812: {
815: da->ops->creatematrix = f;
816: return(0);
817: }
821: /*
822: Creates "balanced" ownership ranges after refinement, constrained by the need for the
823: fine grid boundaries to fall within one stencil width of the coarse partition.
825: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
826: */
827: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
828: {
829: PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
833: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
834: if (ratio == 1) {
835: PetscMemcpy(lf,lc,m*sizeof(lc[0]));
836: return(0);
837: }
838: for (i=0; i<m; i++) totalc += lc[i];
839: remaining = (!periodic) + ratio * (totalc - (!periodic));
840: for (i=0; i<m; i++) {
841: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
842: if (i == m-1) lf[i] = want;
843: else {
844: const PetscInt nextc = startc + lc[i];
845: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
846: * coarse stencil width of the first coarse node in the next subdomain. */
847: while ((startf+want)/ratio < nextc - stencil_width) want++;
848: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
849: * coarse stencil width of the last coarse node in the current subdomain. */
850: while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
851: /* Make sure all constraints are satisfied */
852: if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
853: || ((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");
854: }
855: lf[i] = want;
856: startc += lc[i];
857: startf += lf[i];
858: remaining -= lf[i];
859: }
860: return(0);
861: }
865: /*
866: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
867: fine grid boundaries to fall within one stencil width of the coarse partition.
869: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
870: */
871: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
872: {
873: PetscInt i,totalf,remaining,startc,startf;
877: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
878: if (ratio == 1) {
879: PetscMemcpy(lc,lf,m*sizeof(lf[0]));
880: return(0);
881: }
882: for (i=0,totalf=0; i<m; i++) totalf += lf[i];
883: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
884: for (i=0,startc=0,startf=0; i<m; i++) {
885: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
886: if (i == m-1) lc[i] = want;
887: else {
888: const PetscInt nextf = startf+lf[i];
889: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
890: * node is within one stencil width. */
891: while (nextf/ratio < startc+want-stencil_width) want--;
892: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
893: * fine node is within one stencil width. */
894: while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
895: if (want < 0 || want > remaining
896: || (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");
897: }
898: lc[i] = want;
899: startc += lc[i];
900: startf += lf[i];
901: remaining -= lc[i];
902: }
903: return(0);
904: }
908: PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
909: {
911: PetscInt M,N,P,i;
912: DM da2;
913: DM_DA *dd = (DM_DA*)da->data,*dd2;
919: if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
920: M = dd->refine_x*dd->M;
921: } else {
922: M = 1 + dd->refine_x*(dd->M - 1);
923: }
924: if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
925: if (dd->dim > 1) {
926: N = dd->refine_y*dd->N;
927: } else {
928: N = 1;
929: }
930: } else {
931: N = 1 + dd->refine_y*(dd->N - 1);
932: }
933: if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
934: if (dd->dim > 2) {
935: P = dd->refine_z*dd->P;
936: } else {
937: P = 1;
938: }
939: } else {
940: P = 1 + dd->refine_z*(dd->P - 1);
941: }
942: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
943: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
944: DMDASetDim(da2,dd->dim);
945: DMDASetSizes(da2,M,N,P);
946: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
947: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
948: DMDASetDof(da2,dd->w);
949: DMDASetStencilType(da2,dd->stencil_type);
950: DMDASetStencilWidth(da2,dd->s);
951: if (dd->dim == 3) {
952: PetscInt *lx,*ly,*lz;
953: PetscMalloc3(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly,dd->p,PetscInt,&lz);
954: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
955: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
956: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
957: DMDASetOwnershipRanges(da2,lx,ly,lz);
958: PetscFree3(lx,ly,lz);
959: } else if (dd->dim == 2) {
960: PetscInt *lx,*ly;
961: PetscMalloc2(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly);
962: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
963: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
964: DMDASetOwnershipRanges(da2,lx,ly,NULL);
965: PetscFree2(lx,ly);
966: } else if (dd->dim == 1) {
967: PetscInt *lx;
968: PetscMalloc(dd->m*sizeof(PetscInt),&lx);
969: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
970: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
971: PetscFree(lx);
972: }
973: dd2 = (DM_DA*)da2->data;
975: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
976: da2->ops->creatematrix = da->ops->creatematrix;
977: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
978: da2->ops->getcoloring = da->ops->getcoloring;
979: dd2->interptype = dd->interptype;
981: /* copy fill information if given */
982: if (dd->dfill) {
983: PetscMalloc((dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->dfill);
984: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
985: }
986: if (dd->ofill) {
987: PetscMalloc((dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->ofill);
988: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
989: }
990: /* copy the refine information */
991: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
992: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
993: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
995: /* copy vector type information */
996: PetscFree(da2->vectype);
997: PetscStrallocpy(da->vectype,(char**)&da2->vectype);
999: dd2->lf = dd->lf;
1000: dd2->lj = dd->lj;
1002: da2->leveldown = da->leveldown;
1003: da2->levelup = da->levelup + 1;
1005: DMSetFromOptions(da2);
1006: DMSetUp(da2);
1007: DMViewFromOptions(da2,NULL,"-dm_view");
1009: /* interpolate coordinates if they are set on the coarse grid */
1010: if (da->coordinates) {
1011: DM cdaf,cdac;
1012: Vec coordsc,coordsf;
1013: Mat II;
1015: DMGetCoordinateDM(da,&cdac);
1016: DMGetCoordinates(da,&coordsc);
1017: DMGetCoordinateDM(da2,&cdaf);
1018: /* force creation of the coordinate vector */
1019: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1020: DMGetCoordinates(da2,&coordsf);
1021: DMCreateInterpolation(cdac,cdaf,&II,NULL);
1022: MatInterpolate(II,coordsc,coordsf);
1023: MatDestroy(&II);
1024: }
1026: for (i=0; i<da->bs; i++) {
1027: const char *fieldname;
1028: DMDAGetFieldName(da,i,&fieldname);
1029: DMDASetFieldName(da2,i,fieldname);
1030: }
1032: *daref = da2;
1033: return(0);
1034: }
1039: PetscErrorCode DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1040: {
1042: PetscInt M,N,P,i;
1043: DM da2;
1044: DM_DA *dd = (DM_DA*)da->data,*dd2;
1050: if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1051: M = dd->M / dd->coarsen_x;
1052: } else {
1053: M = 1 + (dd->M - 1) / dd->coarsen_x;
1054: }
1055: if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1056: if (dd->dim > 1) {
1057: N = dd->N / dd->coarsen_y;
1058: } else {
1059: N = 1;
1060: }
1061: } else {
1062: N = 1 + (dd->N - 1) / dd->coarsen_y;
1063: }
1064: if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1065: if (dd->dim > 2) {
1066: P = dd->P / dd->coarsen_z;
1067: } else {
1068: P = 1;
1069: }
1070: } else {
1071: P = 1 + (dd->P - 1) / dd->coarsen_z;
1072: }
1073: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1074: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1075: DMDASetDim(da2,dd->dim);
1076: DMDASetSizes(da2,M,N,P);
1077: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1078: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1079: DMDASetDof(da2,dd->w);
1080: DMDASetStencilType(da2,dd->stencil_type);
1081: DMDASetStencilWidth(da2,dd->s);
1082: if (dd->dim == 3) {
1083: PetscInt *lx,*ly,*lz;
1084: PetscMalloc3(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly,dd->p,PetscInt,&lz);
1085: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1086: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1087: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1088: DMDASetOwnershipRanges(da2,lx,ly,lz);
1089: PetscFree3(lx,ly,lz);
1090: } else if (dd->dim == 2) {
1091: PetscInt *lx,*ly;
1092: PetscMalloc2(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly);
1093: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1094: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1095: DMDASetOwnershipRanges(da2,lx,ly,NULL);
1096: PetscFree2(lx,ly);
1097: } else if (dd->dim == 1) {
1098: PetscInt *lx;
1099: PetscMalloc(dd->m*sizeof(PetscInt),&lx);
1100: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1101: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1102: PetscFree(lx);
1103: }
1104: dd2 = (DM_DA*)da2->data;
1106: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1107: /* da2->ops->createinterpolation = da->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1108: da2->ops->creatematrix = da->ops->creatematrix;
1109: da2->ops->getcoloring = da->ops->getcoloring;
1110: dd2->interptype = dd->interptype;
1112: /* copy fill information if given */
1113: if (dd->dfill) {
1114: PetscMalloc((dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->dfill);
1115: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1116: }
1117: if (dd->ofill) {
1118: PetscMalloc((dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->ofill);
1119: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1120: }
1121: /* copy the refine information */
1122: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1123: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1124: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1126: /* copy vector type information */
1127: PetscFree(da2->vectype);
1128: PetscStrallocpy(da->vectype,(char**)&da2->vectype);
1130: dd2->lf = dd->lf;
1131: dd2->lj = dd->lj;
1133: da2->leveldown = da->leveldown + 1;
1134: da2->levelup = da->levelup;
1136: DMSetFromOptions(da2);
1137: DMSetUp(da2);
1138: DMViewFromOptions(da2,NULL,"-dm_view");
1140: /* inject coordinates if they are set on the fine grid */
1141: if (da->coordinates) {
1142: DM cdaf,cdac;
1143: Vec coordsc,coordsf;
1144: VecScatter inject;
1146: DMGetCoordinateDM(da,&cdaf);
1147: DMGetCoordinates(da,&coordsf);
1148: DMGetCoordinateDM(da2,&cdac);
1149: /* force creation of the coordinate vector */
1150: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1151: DMGetCoordinates(da2,&coordsc);
1153: DMCreateInjection(cdac,cdaf,&inject);
1154: VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1155: VecScatterEnd(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1156: VecScatterDestroy(&inject);
1157: }
1159: for (i=0; i<da->bs; i++) {
1160: const char *fieldname;
1161: DMDAGetFieldName(da,i,&fieldname);
1162: DMDASetFieldName(da2,i,fieldname);
1163: }
1165: *daref = da2;
1166: return(0);
1167: }
1171: PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1172: {
1174: PetscInt i,n,*refx,*refy,*refz;
1178: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1179: if (nlevels == 0) return(0);
1182: /* Get refinement factors, defaults taken from the coarse DMDA */
1183: PetscMalloc3(nlevels,PetscInt,&refx,nlevels,PetscInt,&refy,nlevels,PetscInt,&refz);
1184: for (i=0; i<nlevels; i++) {
1185: DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1186: }
1187: n = nlevels;
1188: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1189: n = nlevels;
1190: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1191: n = nlevels;
1192: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);
1194: DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1195: DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1196: for (i=1; i<nlevels; i++) {
1197: DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1198: DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1199: }
1200: PetscFree3(refx,refy,refz);
1201: return(0);
1202: }
1206: PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1207: {
1209: PetscInt i;
1213: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1214: if (nlevels == 0) return(0);
1216: DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1217: for (i=1; i<nlevels; i++) {
1218: DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1219: }
1220: return(0);
1221: }