Actual source code: da.c
petsc-3.10.5 2019-03-28
1: #include <petsc/private/dmdaimpl.h>
3: /*@
4: DMDASetSizes - Sets the number of grid points in the three dimensional directions
6: Logically Collective on DMDA
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: Since the dimension may not yet have been set the code cannot error check for non-positive Y and Z number of grid points
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()");
30: if (M < 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in X direction must be positive");
31: if (N < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in Y direction must be positive");
32: if (P < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in Z direction must be positive");
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 DMDA
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(), DMDAGetSize(), PetscSplitOwnership()
54: @*/
55: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
56: {
57: DM_DA *dd = (DM_DA*)da->data;
65: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
66: dd->m = m;
67: dd->n = n;
68: dd->p = p;
69: if (da->dim == 2) {
70: PetscMPIInt size;
71: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
72: if ((dd->m > 0) && (dd->n < 0)) {
73: dd->n = size/dd->m;
74: 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);
75: }
76: if ((dd->n > 0) && (dd->m < 0)) {
77: dd->m = size/dd->n;
78: 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);
79: }
80: }
81: return(0);
82: }
84: /*@
85: DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.
87: Not collective
89: Input Parameter:
90: + da - The DMDA
91: - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC
93: Level: intermediate
95: .keywords: distributed array, periodicity
96: .seealso: DMDACreate(), DMDestroy(), DMDA, DMBoundaryType
97: @*/
98: PetscErrorCode DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)
99: {
100: DM_DA *dd = (DM_DA*)da->data;
107: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
108: dd->bx = bx;
109: dd->by = by;
110: dd->bz = bz;
111: return(0);
112: }
114: /*@
115: DMDASetDof - Sets the number of degrees of freedom per vertex
117: Not collective
119: Input Parameters:
120: + da - The DMDA
121: - dof - Number of degrees of freedom
123: Level: intermediate
125: .keywords: distributed array, degrees of freedom
126: .seealso: DMDAGetDof(), DMDACreate(), DMDestroy(), DMDA
127: @*/
128: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
129: {
130: 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->w = dof;
137: da->bs = dof;
138: return(0);
139: }
141: /*@
142: DMDAGetDof - Gets the number of degrees of freedom per vertex
144: Not collective
146: Input Parameter:
147: . da - The DMDA
149: Output Parameter:
150: . dof - Number of degrees of freedom
152: Level: intermediate
154: .keywords: distributed array, degrees of freedom
155: .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
156: @*/
157: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
158: {
159: DM_DA *dd = (DM_DA *) da->data;
164: *dof = dd->w;
165: return(0);
166: }
168: /*@
169: DMDAGetOverlap - Gets the size of the per-processor overlap.
171: Not collective
173: Input Parameters:
174: . da - The DMDA
176: Output Parameters:
177: + x - Overlap in the x direction
178: . y - Overlap in the y direction
179: - z - Overlap in the z direction
181: Level: intermediate
183: .keywords: distributed array, overlap, domain decomposition
184: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
185: @*/
186: PetscErrorCode DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
187: {
188: DM_DA *dd = (DM_DA*)da->data;
192: if (x) *x = dd->xol;
193: if (y) *y = dd->yol;
194: if (z) *z = dd->zol;
195: return(0);
196: }
198: /*@
199: DMDASetOverlap - Sets the size of the per-processor overlap.
201: Not collective
203: Input Parameters:
204: + da - The DMDA
205: . x - Overlap in the x direction
206: . y - Overlap in the y direction
207: - z - Overlap in the z direction
209: Level: intermediate
211: .keywords: distributed array, overlap, domain decomposition
212: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
213: @*/
214: PetscErrorCode DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
215: {
216: DM_DA *dd = (DM_DA*)da->data;
223: dd->xol = x;
224: dd->yol = y;
225: dd->zol = z;
226: return(0);
227: }
230: /*@
231: DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.
233: Not collective
235: Input Parameters:
236: . da - The DMDA
238: Output Parameters:
239: + Nsub - Number of local subdomains created upon decomposition
241: Level: intermediate
243: .keywords: distributed array, domain decomposition
244: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
245: @*/
246: PetscErrorCode DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
247: {
248: DM_DA *dd = (DM_DA*)da->data;
252: if (Nsub) *Nsub = dd->Nsub;
253: return(0);
254: }
256: /*@
257: DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.
259: Not collective
261: Input Parameters:
262: + da - The DMDA
263: - Nsub - The number of local subdomains requested
265: Level: intermediate
267: .keywords: distributed array, domain decomposition
268: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
269: @*/
270: PetscErrorCode DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
271: {
272: DM_DA *dd = (DM_DA*)da->data;
277: dd->Nsub = Nsub;
278: return(0);
279: }
281: /*@
282: DMDASetOffset - Sets the index offset of the DA.
284: Collective on DA
286: Input Parameter:
287: + da - The DMDA
288: . xo - The offset in the x direction
289: . yo - The offset in the y direction
290: - zo - The offset in the z direction
292: Level: intermediate
294: Notes:
295: This is used primarily to overlap a computation on a local DA with that on a global DA without
296: changing boundary conditions or subdomain features that depend upon the global offsets.
298: .keywords: distributed array, degrees of freedom
299: .seealso: DMDAGetOffset(), DMDAVecGetArray()
300: @*/
301: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
302: {
304: DM_DA *dd = (DM_DA*)da->data;
314: dd->xo = xo;
315: dd->yo = yo;
316: dd->zo = zo;
317: dd->Mo = Mo;
318: dd->No = No;
319: dd->Po = Po;
321: if (da->coordinateDM) {
322: DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
323: }
324: return(0);
325: }
327: /*@
328: DMDAGetOffset - Gets the index offset of the DA.
330: Not collective
332: Input Parameter:
333: . da - The DMDA
335: Output Parameters:
336: + xo - The offset in the x direction
337: . yo - The offset in the y direction
338: . zo - The offset in the z direction
339: . Mo - The global size in the x direction
340: . No - The global size in the y direction
341: - Po - The global size in the z direction
343: Level: intermediate
345: .keywords: distributed array, degrees of freedom
346: .seealso: DMDASetOffset(), DMDAVecGetArray()
347: @*/
348: PetscErrorCode DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
349: {
350: DM_DA *dd = (DM_DA*)da->data;
354: if (xo) *xo = dd->xo;
355: if (yo) *yo = dd->yo;
356: if (zo) *zo = dd->zo;
357: if (Mo) *Mo = dd->Mo;
358: if (No) *No = dd->No;
359: if (Po) *Po = dd->Po;
360: return(0);
361: }
363: /*@
364: DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.
366: Not collective
368: Input Parameter:
369: . da - The DMDA
371: Output Parameters:
372: + xs - The start of the region in x
373: . ys - The start of the region in y
374: . zs - The start of the region in z
375: . xs - The size of the region in x
376: . ys - The size of the region in y
377: . zs - The size of the region in z
379: Level: intermediate
381: .keywords: distributed array, degrees of freedom
382: .seealso: DMDAGetOffset(), DMDAVecGetArray()
383: @*/
384: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
385: {
386: DM_DA *dd = (DM_DA*)da->data;
390: if (xs) *xs = dd->nonxs;
391: if (ys) *ys = dd->nonys;
392: if (zs) *zs = dd->nonzs;
393: if (xm) *xm = dd->nonxm;
394: if (ym) *ym = dd->nonym;
395: if (zm) *zm = dd->nonzm;
396: return(0);
397: }
400: /*@
401: DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.
403: Collective on DA
405: Input Parameter:
406: + da - The DMDA
407: . xs - The start of the region in x
408: . ys - The start of the region in y
409: . zs - The start of the region in z
410: . xs - The size of the region in x
411: . ys - The size of the region in y
412: . zs - The size of the region in z
414: Level: intermediate
416: .keywords: distributed array, degrees of freedom
417: .seealso: DMDAGetOffset(), DMDAVecGetArray()
418: @*/
419: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
420: {
421: DM_DA *dd = (DM_DA*)da->data;
431: dd->nonxs = xs;
432: dd->nonys = ys;
433: dd->nonzs = zs;
434: dd->nonxm = xm;
435: dd->nonym = ym;
436: dd->nonzm = zm;
438: return(0);
439: }
441: /*@
442: DMDASetStencilType - Sets the type of the communication stencil
444: Logically Collective on DMDA
446: Input Parameter:
447: + da - The DMDA
448: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
450: Level: intermediate
452: .keywords: distributed array, stencil
453: .seealso: DMDACreate(), DMDestroy(), DMDA
454: @*/
455: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
456: {
457: DM_DA *dd = (DM_DA*)da->data;
462: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
463: dd->stencil_type = stype;
464: return(0);
465: }
467: /*@
468: DMDAGetStencilType - Gets the type of the communication stencil
470: Not collective
472: Input Parameter:
473: . da - The DMDA
475: Output Parameter:
476: . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
478: Level: intermediate
480: .keywords: distributed array, stencil
481: .seealso: DMDACreate(), DMDestroy(), DMDA
482: @*/
483: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
484: {
485: DM_DA *dd = (DM_DA*)da->data;
490: *stype = dd->stencil_type;
491: return(0);
492: }
494: /*@
495: DMDASetStencilWidth - Sets the width of the communication stencil
497: Logically Collective on DMDA
499: Input Parameter:
500: + da - The DMDA
501: - width - The stencil width
503: Level: intermediate
505: .keywords: distributed array, stencil
506: .seealso: DMDACreate(), DMDestroy(), DMDA
507: @*/
508: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
509: {
510: DM_DA *dd = (DM_DA*)da->data;
515: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
516: dd->s = width;
517: return(0);
518: }
520: /*@
521: DMDAGetStencilWidth - Gets the width of the communication stencil
523: Not collective
525: Input Parameter:
526: . da - The DMDA
528: Output Parameter:
529: . width - The stencil width
531: Level: intermediate
533: .keywords: distributed array, stencil
534: .seealso: DMDACreate(), DMDestroy(), DMDA
535: @*/
536: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
537: {
538: DM_DA *dd = (DM_DA *) da->data;
543: *width = dd->s;
544: return(0);
545: }
547: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
548: {
549: PetscInt i,sum;
552: if (M < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
553: for (i=sum=0; i<m; i++) sum += lx[i];
554: if (sum != M) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
555: return(0);
556: }
558: /*@
559: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
561: Logically Collective on DMDA
563: Input Parameter:
564: + da - The DMDA
565: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
566: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
567: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.
569: Level: intermediate
571: Note: these numbers are NOT multiplied by the number of dof per node.
573: .keywords: distributed array
574: .seealso: DMDACreate(), DMDestroy(), DMDA
575: @*/
576: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
577: {
579: DM_DA *dd = (DM_DA*)da->data;
583: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
584: if (lx) {
585: if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
586: DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
587: if (!dd->lx) {
588: PetscMalloc1(dd->m, &dd->lx);
589: }
590: PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
591: }
592: if (ly) {
593: if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
594: DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
595: if (!dd->ly) {
596: PetscMalloc1(dd->n, &dd->ly);
597: }
598: PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
599: }
600: if (lz) {
601: if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
602: DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
603: if (!dd->lz) {
604: PetscMalloc1(dd->p, &dd->lz);
605: }
606: PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
607: }
608: return(0);
609: }
611: /*@
612: DMDASetInterpolationType - Sets the type of interpolation that will be
613: returned by DMCreateInterpolation()
615: Logically Collective on DMDA
617: Input Parameter:
618: + da - initial distributed array
619: . ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
621: Level: intermediate
623: Notes:
624: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
626: .keywords: distributed array, interpolation
628: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
629: @*/
630: PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
631: {
632: DM_DA *dd = (DM_DA*)da->data;
637: dd->interptype = ctype;
638: return(0);
639: }
641: /*@
642: DMDAGetInterpolationType - Gets the type of interpolation that will be
643: used by DMCreateInterpolation()
645: Not Collective
647: Input Parameter:
648: . da - distributed array
650: Output Parameter:
651: . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
653: Level: intermediate
655: .keywords: distributed array, interpolation
657: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
658: @*/
659: PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
660: {
661: DM_DA *dd = (DM_DA*)da->data;
666: *ctype = dd->interptype;
667: return(0);
668: }
670: /*@C
671: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
672: processes neighbors.
674: Not Collective
676: Input Parameter:
677: . da - the DMDA object
679: Output Parameters:
680: . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
681: this process itself is in the list
683: Notes:
684: In 2d the array is of length 9, in 3d of length 27
685: Not supported in 1d
686: Do not free the array, it is freed when the DMDA is destroyed.
688: Fortran Notes:
689: In fortran you must pass in an array of the appropriate length.
691: Level: intermediate
693: @*/
694: PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
695: {
696: DM_DA *dd = (DM_DA*)da->data;
700: *ranks = dd->neighbors;
701: return(0);
702: }
704: /*@C
705: DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
707: Not Collective
709: Input Parameter:
710: . da - the DMDA object
712: Output Parameter:
713: + lx - ownership along x direction (optional)
714: . ly - ownership along y direction (optional)
715: - lz - ownership along z direction (optional)
717: Level: intermediate
719: Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
721: In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
722: eighth arguments from DMDAGetInfo()
724: In C you should not free these arrays, nor change the values in them. They will only have valid values while the
725: DMDA they came from still exists (has not been destroyed).
727: These numbers are NOT multiplied by the number of dof per node.
729: Not available from Fortran
731: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
732: @*/
733: PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
734: {
735: DM_DA *dd = (DM_DA*)da->data;
739: if (lx) *lx = dd->lx;
740: if (ly) *ly = dd->ly;
741: if (lz) *lz = dd->lz;
742: return(0);
743: }
745: /*@
746: DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
748: Logically Collective on DMDA
750: Input Parameters:
751: + da - the DMDA object
752: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
753: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
754: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
756: Options Database:
757: + -da_refine_x - refinement ratio in x direction
758: . -da_refine_y - refinement ratio in y direction
759: - -da_refine_z - refinement ratio in z direction
761: Level: intermediate
763: Notes:
764: Pass PETSC_IGNORE to leave a value unchanged
766: .seealso: DMRefine(), DMDAGetRefinementFactor()
767: @*/
768: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
769: {
770: DM_DA *dd = (DM_DA*)da->data;
778: if (refine_x > 0) dd->refine_x = refine_x;
779: if (refine_y > 0) dd->refine_y = refine_y;
780: if (refine_z > 0) dd->refine_z = refine_z;
781: return(0);
782: }
784: /*@C
785: DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
787: Not Collective
789: Input Parameter:
790: . da - the DMDA object
792: Output Parameters:
793: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
794: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
795: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
797: Level: intermediate
799: Notes:
800: Pass NULL for values you do not need
802: .seealso: DMRefine(), DMDASetRefinementFactor()
803: @*/
804: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
805: {
806: DM_DA *dd = (DM_DA*)da->data;
810: if (refine_x) *refine_x = dd->refine_x;
811: if (refine_y) *refine_y = dd->refine_y;
812: if (refine_z) *refine_z = dd->refine_z;
813: return(0);
814: }
816: /*@C
817: DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
819: Logically Collective on DMDA
821: Input Parameters:
822: + da - the DMDA object
823: - f - the function that allocates the matrix for that specific DMDA
825: Level: developer
827: Notes:
828: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
829: the diagonal and off-diagonal blocks of the matrix
831: Not supported from Fortran
833: .seealso: DMCreateMatrix(), DMDASetBlockFills()
834: @*/
835: PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
836: {
839: da->ops->creatematrix = f;
840: return(0);
841: }
843: /*
844: Creates "balanced" ownership ranges after refinement, constrained by the need for the
845: fine grid boundaries to fall within one stencil width of the coarse partition.
847: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
848: */
849: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
850: {
851: PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
855: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
856: if (ratio == 1) {
857: PetscMemcpy(lf,lc,m*sizeof(lc[0]));
858: return(0);
859: }
860: for (i=0; i<m; i++) totalc += lc[i];
861: remaining = (!periodic) + ratio * (totalc - (!periodic));
862: for (i=0; i<m; i++) {
863: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
864: if (i == m-1) lf[i] = want;
865: else {
866: const PetscInt nextc = startc + lc[i];
867: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
868: * coarse stencil width of the first coarse node in the next subdomain. */
869: while ((startf+want)/ratio < nextc - stencil_width) want++;
870: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
871: * coarse stencil width of the last coarse node in the current subdomain. */
872: while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
873: /* Make sure all constraints are satisfied */
874: if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
875: || ((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");
876: }
877: lf[i] = want;
878: startc += lc[i];
879: startf += lf[i];
880: remaining -= lf[i];
881: }
882: return(0);
883: }
885: /*
886: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
887: fine grid boundaries to fall within one stencil width of the coarse partition.
889: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
890: */
891: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
892: {
893: PetscInt i,totalf,remaining,startc,startf;
897: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
898: if (ratio == 1) {
899: PetscMemcpy(lc,lf,m*sizeof(lf[0]));
900: return(0);
901: }
902: for (i=0,totalf=0; i<m; i++) totalf += lf[i];
903: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
904: for (i=0,startc=0,startf=0; i<m; i++) {
905: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
906: if (i == m-1) lc[i] = want;
907: else {
908: const PetscInt nextf = startf+lf[i];
909: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
910: * node is within one stencil width. */
911: while (nextf/ratio < startc+want-stencil_width) want--;
912: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
913: * fine node is within one stencil width. */
914: while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
915: if (want < 0 || want > remaining
916: || (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");
917: }
918: lc[i] = want;
919: startc += lc[i];
920: startf += lf[i];
921: remaining -= lc[i];
922: }
923: return(0);
924: }
926: PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
927: {
929: PetscInt M,N,P,i,dim;
930: DM da2;
931: DM_DA *dd = (DM_DA*)da->data,*dd2;
937: DMGetDimension(da, &dim);
938: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
939: M = dd->refine_x*dd->M;
940: } else {
941: M = 1 + dd->refine_x*(dd->M - 1);
942: }
943: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
944: if (dim > 1) {
945: N = dd->refine_y*dd->N;
946: } else {
947: N = 1;
948: }
949: } else {
950: N = 1 + dd->refine_y*(dd->N - 1);
951: }
952: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
953: if (dim > 2) {
954: P = dd->refine_z*dd->P;
955: } else {
956: P = 1;
957: }
958: } else {
959: P = 1 + dd->refine_z*(dd->P - 1);
960: }
961: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
962: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
963: DMSetDimension(da2,dim);
964: DMDASetSizes(da2,M,N,P);
965: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
966: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
967: DMDASetDof(da2,dd->w);
968: DMDASetStencilType(da2,dd->stencil_type);
969: DMDASetStencilWidth(da2,dd->s);
970: if (dim == 3) {
971: PetscInt *lx,*ly,*lz;
972: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
973: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
974: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
975: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
976: DMDASetOwnershipRanges(da2,lx,ly,lz);
977: PetscFree3(lx,ly,lz);
978: } else if (dim == 2) {
979: PetscInt *lx,*ly;
980: PetscMalloc2(dd->m,&lx,dd->n,&ly);
981: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
982: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
983: DMDASetOwnershipRanges(da2,lx,ly,NULL);
984: PetscFree2(lx,ly);
985: } else if (dim == 1) {
986: PetscInt *lx;
987: PetscMalloc1(dd->m,&lx);
988: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
989: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
990: PetscFree(lx);
991: }
992: dd2 = (DM_DA*)da2->data;
994: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
995: da2->ops->creatematrix = da->ops->creatematrix;
996: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
997: da2->ops->getcoloring = da->ops->getcoloring;
998: dd2->interptype = dd->interptype;
1000: /* copy fill information if given */
1001: if (dd->dfill) {
1002: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1003: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1004: }
1005: if (dd->ofill) {
1006: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1007: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1008: }
1009: /* copy the refine information */
1010: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
1011: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
1012: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
1014: if (dd->refine_z_hier) {
1015: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) {
1016: dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
1017: }
1018: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) {
1019: dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
1020: }
1021: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1022: PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1023: PetscMemcpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n*sizeof(PetscInt));
1024: }
1025: if (dd->refine_y_hier) {
1026: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) {
1027: dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
1028: }
1029: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) {
1030: dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
1031: }
1032: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1033: PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1034: PetscMemcpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n*sizeof(PetscInt));
1035: }
1036: if (dd->refine_x_hier) {
1037: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) {
1038: dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
1039: }
1040: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) {
1041: dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
1042: }
1043: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1044: PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1045: PetscMemcpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n*sizeof(PetscInt));
1046: }
1049: /* copy vector type information */
1050: DMSetVecType(da2,da->vectype);
1052: dd2->lf = dd->lf;
1053: dd2->lj = dd->lj;
1055: da2->leveldown = da->leveldown;
1056: da2->levelup = da->levelup + 1;
1058: DMSetUp(da2);
1060: /* interpolate coordinates if they are set on the coarse grid */
1061: if (da->coordinates) {
1062: DM cdaf,cdac;
1063: Vec coordsc,coordsf;
1064: Mat II;
1066: DMGetCoordinateDM(da,&cdac);
1067: DMGetCoordinates(da,&coordsc);
1068: DMGetCoordinateDM(da2,&cdaf);
1069: /* force creation of the coordinate vector */
1070: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1071: DMGetCoordinates(da2,&coordsf);
1072: DMCreateInterpolation(cdac,cdaf,&II,NULL);
1073: MatInterpolate(II,coordsc,coordsf);
1074: MatDestroy(&II);
1075: }
1077: for (i=0; i<da->bs; i++) {
1078: const char *fieldname;
1079: DMDAGetFieldName(da,i,&fieldname);
1080: DMDASetFieldName(da2,i,fieldname);
1081: }
1083: *daref = da2;
1084: return(0);
1085: }
1088: PetscErrorCode DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1089: {
1091: PetscInt M,N,P,i,dim;
1092: DM da2;
1093: DM_DA *dd = (DM_DA*)da->data,*dd2;
1099: DMGetDimension(da, &dim);
1100: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1101: M = dd->M / dd->coarsen_x;
1102: } else {
1103: M = 1 + (dd->M - 1) / dd->coarsen_x;
1104: }
1105: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1106: if (dim > 1) {
1107: N = dd->N / dd->coarsen_y;
1108: } else {
1109: N = 1;
1110: }
1111: } else {
1112: N = 1 + (dd->N - 1) / dd->coarsen_y;
1113: }
1114: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1115: if (dim > 2) {
1116: P = dd->P / dd->coarsen_z;
1117: } else {
1118: P = 1;
1119: }
1120: } else {
1121: P = 1 + (dd->P - 1) / dd->coarsen_z;
1122: }
1123: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1124: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1125: DMSetDimension(da2,dim);
1126: DMDASetSizes(da2,M,N,P);
1127: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1128: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1129: DMDASetDof(da2,dd->w);
1130: DMDASetStencilType(da2,dd->stencil_type);
1131: DMDASetStencilWidth(da2,dd->s);
1132: if (dim == 3) {
1133: PetscInt *lx,*ly,*lz;
1134: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1135: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1136: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1137: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1138: DMDASetOwnershipRanges(da2,lx,ly,lz);
1139: PetscFree3(lx,ly,lz);
1140: } else if (dim == 2) {
1141: PetscInt *lx,*ly;
1142: PetscMalloc2(dd->m,&lx,dd->n,&ly);
1143: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1144: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1145: DMDASetOwnershipRanges(da2,lx,ly,NULL);
1146: PetscFree2(lx,ly);
1147: } else if (dim == 1) {
1148: PetscInt *lx;
1149: PetscMalloc1(dd->m,&lx);
1150: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1151: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1152: PetscFree(lx);
1153: }
1154: dd2 = (DM_DA*)da2->data;
1156: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1157: /* da2->ops->createinterpolation = da->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1158: da2->ops->creatematrix = da->ops->creatematrix;
1159: da2->ops->getcoloring = da->ops->getcoloring;
1160: dd2->interptype = dd->interptype;
1162: /* copy fill information if given */
1163: if (dd->dfill) {
1164: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1165: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1166: }
1167: if (dd->ofill) {
1168: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1169: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1170: }
1171: /* copy the refine information */
1172: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1173: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1174: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1176: if (dd->refine_z_hier) {
1177: if (da->levelup - da->leveldown -1 > -1 && da->levelup - da->leveldown - 1< dd->refine_z_hier_n) {
1178: dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown - 1];
1179: }
1180: if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_z_hier_n) {
1181: dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown - 2];
1182: }
1183: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1184: PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1185: PetscMemcpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n*sizeof(PetscInt));
1186: }
1187: if (dd->refine_y_hier) {
1188: if (da->levelup - da->leveldown - 1 > -1 && da->levelup - da->leveldown - 1< dd->refine_y_hier_n) {
1189: dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown - 1];
1190: }
1191: if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_y_hier_n) {
1192: dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown - 2];
1193: }
1194: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1195: PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1196: PetscMemcpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n*sizeof(PetscInt));
1197: }
1198: if (dd->refine_x_hier) {
1199: if (da->levelup - da->leveldown - 1 > -1 && da->levelup - da->leveldown - 1 < dd->refine_x_hier_n) {
1200: dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown - 1];
1201: }
1202: if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_x_hier_n) {
1203: dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown - 2];
1204: }
1205: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1206: PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1207: PetscMemcpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n*sizeof(PetscInt));
1208: }
1210: /* copy vector type information */
1211: DMSetVecType(da2,da->vectype);
1213: dd2->lf = dd->lf;
1214: dd2->lj = dd->lj;
1216: da2->leveldown = da->leveldown + 1;
1217: da2->levelup = da->levelup;
1219: DMSetUp(da2);
1221: /* inject coordinates if they are set on the fine grid */
1222: if (da->coordinates) {
1223: DM cdaf,cdac;
1224: Vec coordsc,coordsf;
1225: Mat inject;
1226: VecScatter vscat;
1228: DMGetCoordinateDM(da,&cdaf);
1229: DMGetCoordinates(da,&coordsf);
1230: DMGetCoordinateDM(da2,&cdac);
1231: /* force creation of the coordinate vector */
1232: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1233: DMGetCoordinates(da2,&coordsc);
1235: DMCreateInjection(cdac,cdaf,&inject);
1236: MatScatterGetVecScatter(inject,&vscat);
1237: VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1238: VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1239: MatDestroy(&inject);
1240: }
1242: for (i=0; i<da->bs; i++) {
1243: const char *fieldname;
1244: DMDAGetFieldName(da,i,&fieldname);
1245: DMDASetFieldName(da2,i,fieldname);
1246: }
1248: *daref = da2;
1249: return(0);
1250: }
1252: PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1253: {
1255: PetscInt i,n,*refx,*refy,*refz;
1259: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1260: if (nlevels == 0) return(0);
1263: /* Get refinement factors, defaults taken from the coarse DMDA */
1264: PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1265: for (i=0; i<nlevels; i++) {
1266: DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1267: }
1268: n = nlevels;
1269: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1270: n = nlevels;
1271: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1272: n = nlevels;
1273: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);
1275: DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1276: DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1277: for (i=1; i<nlevels; i++) {
1278: DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1279: DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1280: }
1281: PetscFree3(refx,refy,refz);
1282: return(0);
1283: }
1285: PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1286: {
1288: PetscInt i;
1292: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1293: if (nlevels == 0) return(0);
1295: DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1296: for (i=1; i<nlevels; i++) {
1297: DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1298: }
1299: return(0);
1300: }
1302: #include <petscgll.h>
1304: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm,PetscGLL *gll)
1305: {
1307: PetscInt i,j,n = gll->n,xs,xn,q;
1308: PetscScalar *xx;
1309: PetscReal h;
1310: Vec x;
1311: DM_DA *da = (DM_DA*)dm->data;
1314: if (da->bx != DM_BOUNDARY_PERIODIC) {
1315: DMDAGetInfo(dm,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1316: q = (q-1)/(n-1); /* number of spectral elements */
1317: h = 2.0/q;
1318: DMDAGetCorners(dm,&xs,NULL,NULL,&xn,NULL,NULL);
1319: xs = xs/(n-1);
1320: xn = xn/(n-1);
1321: DMDASetUniformCoordinates(dm,-1.,1.,0.,0.,0.,0.);
1322: DMGetCoordinates(dm,&x);
1323: DMDAVecGetArray(dm,x,&xx);
1325: /* loop over local spectral elements */
1326: for (j=xs; j<xs+xn; j++) {
1327: /*
1328: Except for the first process, each process starts on the second GLL point of the first element on that process
1329: */
1330: for (i= (j == xs && xs > 0)? 1 : 0; i<n; i++) {
1331: xx[j*(n-1) + i] = -1.0 + h*j + h*(gll->nodes[i]+1.0)/2.;
1332: }
1333: }
1334: DMDAVecRestoreArray(dm,x,&xx);
1335: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for periodic");
1336: return(0);
1337: }
1339: /*@
1341: 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
1343: Collective on DM
1345: Input Parameters:
1346: + da - the DMDA object
1347: - gll - the GLL object
1349: Notes:
1350: the parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1351: on each process much be divisible by the number of GLL elements needed per process. This depends on whether the DM is
1352: periodic or not.
1354: Level: advanced
1356: .seealso: DMDACreate(), PetscGLLCreate(), DMGetCoordinates()
1357: @*/
1358: PetscErrorCode DMDASetGLLCoordinates(DM da,PetscGLL *gll)
1359: {
1363: if (da->dim == 1) {
1364: DMDASetGLLCoordinates_1d(da,gll);
1365: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for 2 or 3d");
1366: return(0);
1367: }
1369: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM da1,DM dm2,PetscBool *compatible,PetscBool *set)
1370: {
1372: DM_DA *dd1 = (DM_DA*)da1->data,*dd2;
1373: DM da2;
1374: DMType dmtype2;
1375: PetscBool isda,compatibleLocal;
1376: PetscInt i;
1379: if (!da1->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da1),PETSC_ERR_ARG_WRONGSTATE,"DMSetUp() must be called on first DM before DMGetCompatibility()");
1380: DMGetType(dm2,&dmtype2);
1381: PetscStrcmp(dmtype2,DMDA,&isda);
1382: if (isda) {
1383: da2 = dm2;
1384: dd2 = (DM_DA*)da2->data;
1385: if (!da2->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da2),PETSC_ERR_ARG_WRONGSTATE,"DMSetUp() must be called on second DM before DMGetCompatibility()");
1386: compatibleLocal = (PetscBool)(da1->dim == da2->dim);
1387: if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
1388: /* Global size ranks Boundary type */
1389: if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
1390: if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
1391: if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
1392: if (compatibleLocal) {
1393: for (i=0; i<dd1->m; ++i) {
1394: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i])); /* Local size */
1395: }
1396: }
1397: if (compatibleLocal && da1->dim > 1) {
1398: for (i=0; i<dd1->n; ++i) {
1399: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
1400: }
1401: }
1402: if (compatibleLocal && da1->dim > 2) {
1403: for (i=0; i<dd1->p; ++i) {
1404: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
1405: }
1406: }
1407: *compatible = compatibleLocal;
1408: *set = PETSC_TRUE;
1409: } else {
1410: /* Decline to determine compatibility with other DM types */
1411: *set = PETSC_FALSE;
1412: }
1413: return(0);
1414: }