Actual source code: da.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
5: /*@
6: DMDASetSizes - Sets the global sizes
8: Logically Collective on DMDA
10: Input Parameters:
11: + da - the DMDA
12: . M - the global X size (or PETSC_DECIDE)
13: . N - the global Y size (or PETSC_DECIDE)
14: - P - the global Z size (or PETSC_DECIDE)
16: Level: intermediate
18: .seealso: DMDAGetSize(), PetscSplitOwnership()
19: @*/
20: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
21: {
22: DM_DA *dd = (DM_DA*)da->data;
29: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
31: dd->M = M;
32: dd->N = N;
33: dd->P = P;
34: return(0);
35: }
39: /*@
40: DMDASetNumProcs - Sets the number of processes in each dimension
42: Logically Collective on DMDA
44: Input Parameters:
45: + da - the DMDA
46: . m - the number of X procs (or PETSC_DECIDE)
47: . n - the number of Y procs (or PETSC_DECIDE)
48: - p - the number of Z procs (or PETSC_DECIDE)
50: Level: intermediate
52: .seealso: DMDASetSizes(), DMDAGetSize(), PetscSplitOwnership()
53: @*/
54: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
55: {
56: DM_DA *dd = (DM_DA*)da->data;
64: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
65: dd->m = m;
66: dd->n = n;
67: dd->p = p;
68: if (da->dim == 2) {
69: PetscMPIInt size;
70: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
71: if ((dd->m > 0) && (dd->n < 0)) {
72: dd->n = size/dd->m;
73: 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);
74: }
75: if ((dd->n > 0) && (dd->m < 0)) {
76: dd->m = size/dd->n;
77: 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);
78: }
79: }
80: return(0);
81: }
85: /*@
86: DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.
88: Not collective
90: Input Parameter:
91: + da - The DMDA
92: - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC
94: Level: intermediate
96: .keywords: distributed array, periodicity
97: .seealso: DMDACreate(), DMDestroy(), DMDA, DMBoundaryType
98: @*/
99: PetscErrorCode DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)
100: {
101: DM_DA *dd = (DM_DA*)da->data;
108: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
109: dd->bx = bx;
110: dd->by = by;
111: dd->bz = bz;
112: return(0);
113: }
117: /*@
118: DMDASetDof - Sets the number of degrees of freedom per vertex
120: Not collective
122: Input Parameters:
123: + da - The DMDA
124: - dof - Number of degrees of freedom
126: Level: intermediate
128: .keywords: distributed array, degrees of freedom
129: .seealso: DMDAGetDof(), DMDACreate(), DMDestroy(), DMDA
130: @*/
131: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
132: {
133: DM_DA *dd = (DM_DA*)da->data;
138: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
139: dd->w = dof;
140: da->bs = dof;
141: return(0);
142: }
146: /*@
147: DMDAGetDof - Gets the number of degrees of freedom per vertex
149: Not collective
151: Input Parameter:
152: . da - The DMDA
154: Output Parameter:
155: . dof - Number of degrees of freedom
157: Level: intermediate
159: .keywords: distributed array, degrees of freedom
160: .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
161: @*/
162: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
163: {
164: DM_DA *dd = (DM_DA *) da->data;
169: *dof = dd->w;
170: return(0);
171: }
175: /*@
176: DMDAGetOverlap - Gets the size of the per-processor overlap.
178: Not collective
180: Input Parameters:
181: . da - The DMDA
183: Output Parameters:
184: + x - Overlap in the x direction
185: . y - Overlap in the y direction
186: - z - Overlap in the z direction
188: Level: intermediate
190: .keywords: distributed array, overlap, domain decomposition
191: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
192: @*/
193: PetscErrorCode DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
194: {
195: DM_DA *dd = (DM_DA*)da->data;
199: if (x) *x = dd->xol;
200: if (y) *y = dd->yol;
201: if (z) *z = dd->zol;
202: return(0);
203: }
207: /*@
208: DMDASetOverlap - Sets the size of the per-processor overlap.
210: Not collective
212: Input Parameters:
213: + da - The DMDA
214: . x - Overlap in the x direction
215: . y - Overlap in the y direction
216: - z - Overlap in the z direction
218: Level: intermediate
220: .keywords: distributed array, overlap, domain decomposition
221: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
222: @*/
223: PetscErrorCode DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
224: {
225: DM_DA *dd = (DM_DA*)da->data;
232: dd->xol = x;
233: dd->yol = y;
234: dd->zol = z;
235: return(0);
236: }
241: /*@
242: DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.
244: Not collective
246: Input Parameters:
247: . da - The DMDA
249: Output Parameters:
250: + Nsub - Number of local subdomains created upon decomposition
252: Level: intermediate
254: .keywords: distributed array, domain decomposition
255: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
256: @*/
257: PetscErrorCode DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
258: {
259: DM_DA *dd = (DM_DA*)da->data;
263: if (Nsub) *Nsub = dd->Nsub;
264: return(0);
265: }
269: /*@
270: DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.
272: Not collective
274: Input Parameters:
275: + da - The DMDA
276: - Nsub - The number of local subdomains requested
278: Level: intermediate
280: .keywords: distributed array, domain decomposition
281: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
282: @*/
283: PetscErrorCode DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
284: {
285: DM_DA *dd = (DM_DA*)da->data;
290: dd->Nsub = Nsub;
291: return(0);
292: }
296: /*@
297: DMDASetOffset - Sets the index offset of the DA.
299: Collective on DA
301: Input Parameter:
302: + da - The DMDA
303: . xo - The offset in the x direction
304: . yo - The offset in the y direction
305: - zo - The offset in the z direction
307: Level: intermediate
309: Notes: This is used primarily to overlap a computation on a local DA with that on a global DA without
310: changing boundary conditions or subdomain features that depend upon the global offsets.
312: .keywords: distributed array, degrees of freedom
313: .seealso: DMDAGetOffset(), DMDAVecGetArray()
314: @*/
315: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
316: {
318: DM_DA *dd = (DM_DA*)da->data;
328: dd->xo = xo;
329: dd->yo = yo;
330: dd->zo = zo;
331: dd->Mo = Mo;
332: dd->No = No;
333: dd->Po = Po;
335: if (da->coordinateDM) {
336: DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
337: }
338: return(0);
339: }
343: /*@
344: DMDAGetOffset - Gets the index offset of the DA.
346: Not collective
348: Input Parameter:
349: . da - The DMDA
351: Output Parameters:
352: + xo - The offset in the x direction
353: . yo - The offset in the y direction
354: . zo - The offset in the z direction
355: . Mo - The global size in the x direction
356: . No - The global size in the y direction
357: - Po - The global size in the z direction
359: Level: intermediate
361: .keywords: distributed array, degrees of freedom
362: .seealso: DMDASetOffset(), DMDAVecGetArray()
363: @*/
364: PetscErrorCode DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
365: {
366: DM_DA *dd = (DM_DA*)da->data;
370: if (xo) *xo = dd->xo;
371: if (yo) *yo = dd->yo;
372: if (zo) *zo = dd->zo;
373: if (Mo) *Mo = dd->Mo;
374: if (No) *No = dd->No;
375: if (Po) *Po = dd->Po;
376: return(0);
377: }
381: /*@
382: DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.
384: Not collective
386: Input Parameter:
387: . da - The DMDA
389: Output Parameters:
390: + xs - The start of the region in x
391: . ys - The start of the region in y
392: . zs - The start of the region in z
393: . xs - The size of the region in x
394: . ys - The size of the region in y
395: . zs - The size of the region in z
397: Level: intermediate
399: .keywords: distributed array, degrees of freedom
400: .seealso: DMDAGetOffset(), DMDAVecGetArray()
401: @*/
402: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
403: {
404: DM_DA *dd = (DM_DA*)da->data;
408: if (xs) *xs = dd->nonxs;
409: if (ys) *ys = dd->nonys;
410: if (zs) *zs = dd->nonzs;
411: if (xm) *xm = dd->nonxm;
412: if (ym) *ym = dd->nonym;
413: if (zm) *zm = dd->nonzm;
414: return(0);
415: }
420: /*@
421: DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.
423: Collective on DA
425: Input Parameter:
426: + da - The DMDA
427: . xs - The start of the region in x
428: . ys - The start of the region in y
429: . zs - The start of the region in z
430: . xs - The size of the region in x
431: . ys - The size of the region in y
432: . zs - The size of the region in z
434: Level: intermediate
436: .keywords: distributed array, degrees of freedom
437: .seealso: DMDAGetOffset(), DMDAVecGetArray()
438: @*/
439: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
440: {
441: DM_DA *dd = (DM_DA*)da->data;
451: dd->nonxs = xs;
452: dd->nonys = ys;
453: dd->nonzs = zs;
454: dd->nonxm = xm;
455: dd->nonym = ym;
456: dd->nonzm = zm;
458: return(0);
459: }
463: /*@
464: DMDASetStencilType - Sets the type of the communication stencil
466: Logically Collective on DMDA
468: Input Parameter:
469: + da - The DMDA
470: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
472: Level: intermediate
474: .keywords: distributed array, stencil
475: .seealso: DMDACreate(), DMDestroy(), DMDA
476: @*/
477: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
478: {
479: DM_DA *dd = (DM_DA*)da->data;
484: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
485: dd->stencil_type = stype;
486: return(0);
487: }
491: /*@
492: DMDAGetStencilType - Gets the type of the communication stencil
494: Not collective
496: Input Parameter:
497: . da - The DMDA
499: Output Parameter:
500: . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
502: Level: intermediate
504: .keywords: distributed array, stencil
505: .seealso: DMDACreate(), DMDestroy(), DMDA
506: @*/
507: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
508: {
509: DM_DA *dd = (DM_DA*)da->data;
514: *stype = dd->stencil_type;
515: return(0);
516: }
520: /*@
521: DMDASetStencilWidth - Sets the width of the communication stencil
523: Logically Collective on DMDA
525: Input Parameter:
526: + da - The DMDA
527: - width - The stencil width
529: Level: intermediate
531: .keywords: distributed array, stencil
532: .seealso: DMDACreate(), DMDestroy(), DMDA
533: @*/
534: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
535: {
536: DM_DA *dd = (DM_DA*)da->data;
541: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
542: dd->s = width;
543: return(0);
544: }
548: /*@
549: DMDAGetStencilWidth - Gets the width of the communication stencil
551: Not collective
553: Input Parameter:
554: . da - The DMDA
556: Output Parameter:
557: . width - The stencil width
559: Level: intermediate
561: .keywords: distributed array, stencil
562: .seealso: DMDACreate(), DMDestroy(), DMDA
563: @*/
564: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
565: {
566: DM_DA *dd = (DM_DA *) da->data;
571: *width = dd->s;
572: return(0);
573: }
577: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
578: {
579: PetscInt i,sum;
582: if (M < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
583: for (i=sum=0; i<m; i++) sum += lx[i];
584: if (sum != M) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
585: return(0);
586: }
590: /*@
591: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
593: Logically Collective on DMDA
595: Input Parameter:
596: + da - The DMDA
597: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
598: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
599: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.
601: Level: intermediate
603: Note: these numbers are NOT multiplied by the number of dof per node.
605: .keywords: distributed array
606: .seealso: DMDACreate(), DMDestroy(), DMDA
607: @*/
608: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
609: {
611: DM_DA *dd = (DM_DA*)da->data;
615: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
616: if (lx) {
617: if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
618: DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
619: if (!dd->lx) {
620: PetscMalloc1(dd->m, &dd->lx);
621: }
622: PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
623: }
624: if (ly) {
625: if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
626: DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
627: if (!dd->ly) {
628: PetscMalloc1(dd->n, &dd->ly);
629: }
630: PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
631: }
632: if (lz) {
633: if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
634: DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
635: if (!dd->lz) {
636: PetscMalloc1(dd->p, &dd->lz);
637: }
638: PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
639: }
640: return(0);
641: }
645: /*@
646: DMDASetInterpolationType - Sets the type of interpolation that will be
647: returned by DMCreateInterpolation()
649: Logically Collective on DMDA
651: Input Parameter:
652: + da - initial distributed array
653: . ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
655: Level: intermediate
657: Notes: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
659: .keywords: distributed array, interpolation
661: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
662: @*/
663: PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
664: {
665: DM_DA *dd = (DM_DA*)da->data;
670: dd->interptype = ctype;
671: return(0);
672: }
676: /*@
677: DMDAGetInterpolationType - Gets the type of interpolation that will be
678: used by DMCreateInterpolation()
680: Not Collective
682: Input Parameter:
683: . da - distributed array
685: Output Parameter:
686: . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
688: Level: intermediate
690: .keywords: distributed array, interpolation
692: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
693: @*/
694: PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
695: {
696: DM_DA *dd = (DM_DA*)da->data;
701: *ctype = dd->interptype;
702: return(0);
703: }
707: /*@C
708: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
709: processes neighbors.
711: Not Collective
713: Input Parameter:
714: . da - the DMDA object
716: Output Parameters:
717: . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
718: this process itself is in the list
720: Notes: In 2d the array is of length 9, in 3d of length 27
721: Not supported in 1d
722: Do not free the array, it is freed when the DMDA is destroyed.
724: Fortran Notes: In fortran you must pass in an array of the appropriate length.
726: Level: intermediate
728: @*/
729: PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
730: {
731: DM_DA *dd = (DM_DA*)da->data;
735: *ranks = dd->neighbors;
736: return(0);
737: }
741: /*@C
742: DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
744: Not Collective
746: Input Parameter:
747: . da - the DMDA object
749: Output Parameter:
750: + lx - ownership along x direction (optional)
751: . ly - ownership along y direction (optional)
752: - lz - ownership along z direction (optional)
754: Level: intermediate
756: Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
758: In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
759: eighth arguments from DMDAGetInfo()
761: In C you should not free these arrays, nor change the values in them. They will only have valid values while the
762: DMDA they came from still exists (has not been destroyed).
764: These numbers are NOT multiplied by the number of dof per node.
766: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
767: @*/
768: PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
769: {
770: DM_DA *dd = (DM_DA*)da->data;
774: if (lx) *lx = dd->lx;
775: if (ly) *ly = dd->ly;
776: if (lz) *lz = dd->lz;
777: return(0);
778: }
782: /*@
783: DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
785: Logically Collective on DMDA
787: Input Parameters:
788: + da - the DMDA object
789: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
790: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
791: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
793: Options Database:
794: + -da_refine_x - refinement ratio in x direction
795: . -da_refine_y - refinement ratio in y direction
796: - -da_refine_z - refinement ratio in z direction
798: Level: intermediate
800: Notes: Pass PETSC_IGNORE to leave a value unchanged
802: .seealso: DMRefine(), DMDAGetRefinementFactor()
803: @*/
804: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
805: {
806: DM_DA *dd = (DM_DA*)da->data;
814: if (refine_x > 0) dd->refine_x = refine_x;
815: if (refine_y > 0) dd->refine_y = refine_y;
816: if (refine_z > 0) dd->refine_z = refine_z;
817: return(0);
818: }
822: /*@C
823: DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
825: Not Collective
827: Input Parameter:
828: . da - the DMDA object
830: Output Parameters:
831: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
832: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
833: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
835: Level: intermediate
837: Notes: Pass NULL for values you do not need
839: .seealso: DMRefine(), DMDASetRefinementFactor()
840: @*/
841: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
842: {
843: DM_DA *dd = (DM_DA*)da->data;
847: if (refine_x) *refine_x = dd->refine_x;
848: if (refine_y) *refine_y = dd->refine_y;
849: if (refine_z) *refine_z = dd->refine_z;
850: return(0);
851: }
855: /*@C
856: DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
858: Logically Collective on DMDA
860: Input Parameters:
861: + da - the DMDA object
862: - f - the function that allocates the matrix for that specific DMDA
864: Level: developer
866: Notes: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
867: the diagonal and off-diagonal blocks of the matrix
869: .seealso: DMCreateMatrix(), DMDASetBlockFills()
870: @*/
871: PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
872: {
875: da->ops->creatematrix = f;
876: return(0);
877: }
881: /*
882: Creates "balanced" ownership ranges after refinement, constrained by the need for the
883: fine grid boundaries to fall within one stencil width of the coarse partition.
885: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
886: */
887: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
888: {
889: PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
893: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
894: if (ratio == 1) {
895: PetscMemcpy(lf,lc,m*sizeof(lc[0]));
896: return(0);
897: }
898: for (i=0; i<m; i++) totalc += lc[i];
899: remaining = (!periodic) + ratio * (totalc - (!periodic));
900: for (i=0; i<m; i++) {
901: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
902: if (i == m-1) lf[i] = want;
903: else {
904: const PetscInt nextc = startc + lc[i];
905: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
906: * coarse stencil width of the first coarse node in the next subdomain. */
907: while ((startf+want)/ratio < nextc - stencil_width) want++;
908: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
909: * coarse stencil width of the last coarse node in the current subdomain. */
910: while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
911: /* Make sure all constraints are satisfied */
912: if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
913: || ((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");
914: }
915: lf[i] = want;
916: startc += lc[i];
917: startf += lf[i];
918: remaining -= lf[i];
919: }
920: return(0);
921: }
925: /*
926: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
927: fine grid boundaries to fall within one stencil width of the coarse partition.
929: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
930: */
931: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
932: {
933: PetscInt i,totalf,remaining,startc,startf;
937: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
938: if (ratio == 1) {
939: PetscMemcpy(lc,lf,m*sizeof(lf[0]));
940: return(0);
941: }
942: for (i=0,totalf=0; i<m; i++) totalf += lf[i];
943: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
944: for (i=0,startc=0,startf=0; i<m; i++) {
945: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
946: if (i == m-1) lc[i] = want;
947: else {
948: const PetscInt nextf = startf+lf[i];
949: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
950: * node is within one stencil width. */
951: while (nextf/ratio < startc+want-stencil_width) want--;
952: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
953: * fine node is within one stencil width. */
954: while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
955: if (want < 0 || want > remaining
956: || (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");
957: }
958: lc[i] = want;
959: startc += lc[i];
960: startf += lf[i];
961: remaining -= lc[i];
962: }
963: return(0);
964: }
968: PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
969: {
971: PetscInt M,N,P,i,dim;
972: DM da2;
973: DM_DA *dd = (DM_DA*)da->data,*dd2;
979: DMGetDimension(da, &dim);
980: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
981: M = dd->refine_x*dd->M;
982: } else {
983: M = 1 + dd->refine_x*(dd->M - 1);
984: }
985: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
986: if (dim > 1) {
987: N = dd->refine_y*dd->N;
988: } else {
989: N = 1;
990: }
991: } else {
992: N = 1 + dd->refine_y*(dd->N - 1);
993: }
994: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
995: if (dim > 2) {
996: P = dd->refine_z*dd->P;
997: } else {
998: P = 1;
999: }
1000: } else {
1001: P = 1 + dd->refine_z*(dd->P - 1);
1002: }
1003: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1004: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1005: DMSetDimension(da2,dim);
1006: DMDASetSizes(da2,M,N,P);
1007: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1008: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1009: DMDASetDof(da2,dd->w);
1010: DMDASetStencilType(da2,dd->stencil_type);
1011: DMDASetStencilWidth(da2,dd->s);
1012: if (dim == 3) {
1013: PetscInt *lx,*ly,*lz;
1014: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1015: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
1016: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
1017: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
1018: DMDASetOwnershipRanges(da2,lx,ly,lz);
1019: PetscFree3(lx,ly,lz);
1020: } else if (dim == 2) {
1021: PetscInt *lx,*ly;
1022: PetscMalloc2(dd->m,&lx,dd->n,&ly);
1023: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
1024: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
1025: DMDASetOwnershipRanges(da2,lx,ly,NULL);
1026: PetscFree2(lx,ly);
1027: } else if (dim == 1) {
1028: PetscInt *lx;
1029: PetscMalloc1(dd->m,&lx);
1030: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
1031: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1032: PetscFree(lx);
1033: }
1034: dd2 = (DM_DA*)da2->data;
1036: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
1037: da2->ops->creatematrix = da->ops->creatematrix;
1038: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
1039: da2->ops->getcoloring = da->ops->getcoloring;
1040: dd2->interptype = dd->interptype;
1042: /* copy fill information if given */
1043: if (dd->dfill) {
1044: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1045: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1046: }
1047: if (dd->ofill) {
1048: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1049: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1050: }
1051: /* copy the refine information */
1052: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
1053: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
1054: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
1056: /* copy vector type information */
1057: DMSetVecType(da2,da->vectype);
1059: dd2->lf = dd->lf;
1060: dd2->lj = dd->lj;
1062: da2->leveldown = da->leveldown;
1063: da2->levelup = da->levelup + 1;
1065: DMSetFromOptions(da2);
1066: DMSetUp(da2);
1068: /* interpolate coordinates if they are set on the coarse grid */
1069: if (da->coordinates) {
1070: DM cdaf,cdac;
1071: Vec coordsc,coordsf;
1072: Mat II;
1074: DMGetCoordinateDM(da,&cdac);
1075: DMGetCoordinates(da,&coordsc);
1076: DMGetCoordinateDM(da2,&cdaf);
1077: /* force creation of the coordinate vector */
1078: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1079: DMGetCoordinates(da2,&coordsf);
1080: DMCreateInterpolation(cdac,cdaf,&II,NULL);
1081: MatInterpolate(II,coordsc,coordsf);
1082: MatDestroy(&II);
1083: }
1085: for (i=0; i<da->bs; i++) {
1086: const char *fieldname;
1087: DMDAGetFieldName(da,i,&fieldname);
1088: DMDASetFieldName(da2,i,fieldname);
1089: }
1091: *daref = da2;
1092: return(0);
1093: }
1098: PetscErrorCode DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1099: {
1101: PetscInt M,N,P,i,dim;
1102: DM da2;
1103: DM_DA *dd = (DM_DA*)da->data,*dd2;
1109: DMGetDimension(da, &dim);
1110: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1111: M = dd->M / dd->coarsen_x;
1112: } else {
1113: M = 1 + (dd->M - 1) / dd->coarsen_x;
1114: }
1115: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1116: if (dim > 1) {
1117: N = dd->N / dd->coarsen_y;
1118: } else {
1119: N = 1;
1120: }
1121: } else {
1122: N = 1 + (dd->N - 1) / dd->coarsen_y;
1123: }
1124: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1125: if (dim > 2) {
1126: P = dd->P / dd->coarsen_z;
1127: } else {
1128: P = 1;
1129: }
1130: } else {
1131: P = 1 + (dd->P - 1) / dd->coarsen_z;
1132: }
1133: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1134: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1135: DMSetDimension(da2,dim);
1136: DMDASetSizes(da2,M,N,P);
1137: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1138: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1139: DMDASetDof(da2,dd->w);
1140: DMDASetStencilType(da2,dd->stencil_type);
1141: DMDASetStencilWidth(da2,dd->s);
1142: if (dim == 3) {
1143: PetscInt *lx,*ly,*lz;
1144: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1145: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1146: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1147: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1148: DMDASetOwnershipRanges(da2,lx,ly,lz);
1149: PetscFree3(lx,ly,lz);
1150: } else if (dim == 2) {
1151: PetscInt *lx,*ly;
1152: PetscMalloc2(dd->m,&lx,dd->n,&ly);
1153: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1154: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1155: DMDASetOwnershipRanges(da2,lx,ly,NULL);
1156: PetscFree2(lx,ly);
1157: } else if (dim == 1) {
1158: PetscInt *lx;
1159: PetscMalloc1(dd->m,&lx);
1160: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1161: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1162: PetscFree(lx);
1163: }
1164: dd2 = (DM_DA*)da2->data;
1166: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1167: /* da2->ops->createinterpolation = da->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1168: da2->ops->creatematrix = da->ops->creatematrix;
1169: da2->ops->getcoloring = da->ops->getcoloring;
1170: dd2->interptype = dd->interptype;
1172: /* copy fill information if given */
1173: if (dd->dfill) {
1174: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1175: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1176: }
1177: if (dd->ofill) {
1178: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1179: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1180: }
1181: /* copy the refine information */
1182: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1183: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1184: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1186: /* copy vector type information */
1187: DMSetVecType(da2,da->vectype);
1189: dd2->lf = dd->lf;
1190: dd2->lj = dd->lj;
1192: da2->leveldown = da->leveldown + 1;
1193: da2->levelup = da->levelup;
1195: DMSetFromOptions(da2);
1196: DMSetUp(da2);
1198: /* inject coordinates if they are set on the fine grid */
1199: if (da->coordinates) {
1200: DM cdaf,cdac;
1201: Vec coordsc,coordsf;
1202: Mat inject;
1203: VecScatter vscat;
1205: DMGetCoordinateDM(da,&cdaf);
1206: DMGetCoordinates(da,&coordsf);
1207: DMGetCoordinateDM(da2,&cdac);
1208: /* force creation of the coordinate vector */
1209: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1210: DMGetCoordinates(da2,&coordsc);
1212: DMCreateInjection(cdac,cdaf,&inject);
1213: MatScatterGetVecScatter(inject,&vscat);
1214: VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1215: VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1216: MatDestroy(&inject);
1217: }
1219: for (i=0; i<da->bs; i++) {
1220: const char *fieldname;
1221: DMDAGetFieldName(da,i,&fieldname);
1222: DMDASetFieldName(da2,i,fieldname);
1223: }
1225: *daref = da2;
1226: return(0);
1227: }
1231: PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1232: {
1234: PetscInt i,n,*refx,*refy,*refz;
1238: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1239: if (nlevels == 0) return(0);
1242: /* Get refinement factors, defaults taken from the coarse DMDA */
1243: PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1244: for (i=0; i<nlevels; i++) {
1245: DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1246: }
1247: n = nlevels;
1248: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1249: n = nlevels;
1250: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1251: n = nlevels;
1252: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);
1254: DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1255: DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1256: for (i=1; i<nlevels; i++) {
1257: DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1258: DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1259: }
1260: PetscFree3(refx,refy,refz);
1261: return(0);
1262: }
1266: PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1267: {
1269: PetscInt i;
1273: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1274: if (nlevels == 0) return(0);
1276: DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1277: for (i=1; i<nlevels; i++) {
1278: DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1279: }
1280: return(0);
1281: }