Actual source code: da.c
petsc-3.9.4 2018-09-11
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: This is used primarily to overlap a computation on a local DA with that on a global DA without
295: changing boundary conditions or subdomain features that depend upon the global offsets.
297: .keywords: distributed array, degrees of freedom
298: .seealso: DMDAGetOffset(), DMDAVecGetArray()
299: @*/
300: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
301: {
303: DM_DA *dd = (DM_DA*)da->data;
313: dd->xo = xo;
314: dd->yo = yo;
315: dd->zo = zo;
316: dd->Mo = Mo;
317: dd->No = No;
318: dd->Po = Po;
320: if (da->coordinateDM) {
321: DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
322: }
323: return(0);
324: }
326: /*@
327: DMDAGetOffset - Gets the index offset of the DA.
329: Not collective
331: Input Parameter:
332: . da - The DMDA
334: Output Parameters:
335: + xo - The offset in the x direction
336: . yo - The offset in the y direction
337: . zo - The offset in the z direction
338: . Mo - The global size in the x direction
339: . No - The global size in the y direction
340: - Po - The global size in the z direction
342: Level: intermediate
344: .keywords: distributed array, degrees of freedom
345: .seealso: DMDASetOffset(), DMDAVecGetArray()
346: @*/
347: PetscErrorCode DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
348: {
349: DM_DA *dd = (DM_DA*)da->data;
353: if (xo) *xo = dd->xo;
354: if (yo) *yo = dd->yo;
355: if (zo) *zo = dd->zo;
356: if (Mo) *Mo = dd->Mo;
357: if (No) *No = dd->No;
358: if (Po) *Po = dd->Po;
359: return(0);
360: }
362: /*@
363: DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.
365: Not collective
367: Input Parameter:
368: . da - The DMDA
370: Output Parameters:
371: + xs - The start of the region in x
372: . ys - The start of the region in y
373: . zs - The start of the region in z
374: . xs - The size of the region in x
375: . ys - The size of the region in y
376: . zs - The size of the region in z
378: Level: intermediate
380: .keywords: distributed array, degrees of freedom
381: .seealso: DMDAGetOffset(), DMDAVecGetArray()
382: @*/
383: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
384: {
385: DM_DA *dd = (DM_DA*)da->data;
389: if (xs) *xs = dd->nonxs;
390: if (ys) *ys = dd->nonys;
391: if (zs) *zs = dd->nonzs;
392: if (xm) *xm = dd->nonxm;
393: if (ym) *ym = dd->nonym;
394: if (zm) *zm = dd->nonzm;
395: return(0);
396: }
399: /*@
400: DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.
402: Collective on DA
404: Input Parameter:
405: + da - The DMDA
406: . xs - The start of the region in x
407: . ys - The start of the region in y
408: . zs - The start of the region in z
409: . xs - The size of the region in x
410: . ys - The size of the region in y
411: . zs - The size of the region in z
413: Level: intermediate
415: .keywords: distributed array, degrees of freedom
416: .seealso: DMDAGetOffset(), DMDAVecGetArray()
417: @*/
418: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
419: {
420: DM_DA *dd = (DM_DA*)da->data;
430: dd->nonxs = xs;
431: dd->nonys = ys;
432: dd->nonzs = zs;
433: dd->nonxm = xm;
434: dd->nonym = ym;
435: dd->nonzm = zm;
437: return(0);
438: }
440: /*@
441: DMDASetStencilType - Sets the type of the communication stencil
443: Logically Collective on DMDA
445: Input Parameter:
446: + da - The DMDA
447: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
449: Level: intermediate
451: .keywords: distributed array, stencil
452: .seealso: DMDACreate(), DMDestroy(), DMDA
453: @*/
454: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
455: {
456: DM_DA *dd = (DM_DA*)da->data;
461: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
462: dd->stencil_type = stype;
463: return(0);
464: }
466: /*@
467: DMDAGetStencilType - Gets the type of the communication stencil
469: Not collective
471: Input Parameter:
472: . da - The DMDA
474: Output Parameter:
475: . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
477: Level: intermediate
479: .keywords: distributed array, stencil
480: .seealso: DMDACreate(), DMDestroy(), DMDA
481: @*/
482: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
483: {
484: DM_DA *dd = (DM_DA*)da->data;
489: *stype = dd->stencil_type;
490: return(0);
491: }
493: /*@
494: DMDASetStencilWidth - Sets the width of the communication stencil
496: Logically Collective on DMDA
498: Input Parameter:
499: + da - The DMDA
500: - width - The stencil width
502: Level: intermediate
504: .keywords: distributed array, stencil
505: .seealso: DMDACreate(), DMDestroy(), DMDA
506: @*/
507: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
508: {
509: DM_DA *dd = (DM_DA*)da->data;
514: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
515: dd->s = width;
516: return(0);
517: }
519: /*@
520: DMDAGetStencilWidth - Gets the width of the communication stencil
522: Not collective
524: Input Parameter:
525: . da - The DMDA
527: Output Parameter:
528: . width - The stencil width
530: Level: intermediate
532: .keywords: distributed array, stencil
533: .seealso: DMDACreate(), DMDestroy(), DMDA
534: @*/
535: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
536: {
537: DM_DA *dd = (DM_DA *) da->data;
542: *width = dd->s;
543: return(0);
544: }
546: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
547: {
548: PetscInt i,sum;
551: if (M < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
552: for (i=sum=0; i<m; i++) sum += lx[i];
553: if (sum != M) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
554: return(0);
555: }
557: /*@
558: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
560: Logically Collective on DMDA
562: Input Parameter:
563: + da - The DMDA
564: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
565: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
566: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.
568: Level: intermediate
570: Note: these numbers are NOT multiplied by the number of dof per node.
572: .keywords: distributed array
573: .seealso: DMDACreate(), DMDestroy(), DMDA
574: @*/
575: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
576: {
578: DM_DA *dd = (DM_DA*)da->data;
582: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
583: if (lx) {
584: if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
585: DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
586: if (!dd->lx) {
587: PetscMalloc1(dd->m, &dd->lx);
588: }
589: PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
590: }
591: if (ly) {
592: if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
593: DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
594: if (!dd->ly) {
595: PetscMalloc1(dd->n, &dd->ly);
596: }
597: PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
598: }
599: if (lz) {
600: if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
601: DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
602: if (!dd->lz) {
603: PetscMalloc1(dd->p, &dd->lz);
604: }
605: PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
606: }
607: return(0);
608: }
610: /*@
611: DMDASetInterpolationType - Sets the type of interpolation that will be
612: returned by DMCreateInterpolation()
614: Logically Collective on DMDA
616: Input Parameter:
617: + da - initial distributed array
618: . ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
620: Level: intermediate
622: Notes: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
624: .keywords: distributed array, interpolation
626: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
627: @*/
628: PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
629: {
630: DM_DA *dd = (DM_DA*)da->data;
635: dd->interptype = ctype;
636: return(0);
637: }
639: /*@
640: DMDAGetInterpolationType - Gets the type of interpolation that will be
641: used by DMCreateInterpolation()
643: Not Collective
645: Input Parameter:
646: . da - distributed array
648: Output Parameter:
649: . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
651: Level: intermediate
653: .keywords: distributed array, interpolation
655: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
656: @*/
657: PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
658: {
659: DM_DA *dd = (DM_DA*)da->data;
664: *ctype = dd->interptype;
665: return(0);
666: }
668: /*@C
669: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
670: processes neighbors.
672: Not Collective
674: Input Parameter:
675: . da - the DMDA object
677: Output Parameters:
678: . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
679: this process itself is in the list
681: Notes: In 2d the array is of length 9, in 3d of length 27
682: Not supported in 1d
683: Do not free the array, it is freed when the DMDA is destroyed.
685: Fortran Notes: In fortran you must pass in an array of the appropriate length.
687: Level: intermediate
689: @*/
690: PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
691: {
692: DM_DA *dd = (DM_DA*)da->data;
696: *ranks = dd->neighbors;
697: return(0);
698: }
700: /*@C
701: DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
703: Not Collective
705: Input Parameter:
706: . da - the DMDA object
708: Output Parameter:
709: + lx - ownership along x direction (optional)
710: . ly - ownership along y direction (optional)
711: - lz - ownership along z direction (optional)
713: Level: intermediate
715: Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
717: In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
718: eighth arguments from DMDAGetInfo()
720: In C you should not free these arrays, nor change the values in them. They will only have valid values while the
721: DMDA they came from still exists (has not been destroyed).
723: These numbers are NOT multiplied by the number of dof per node.
725: Not available from Fortran
727: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
728: @*/
729: PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
730: {
731: DM_DA *dd = (DM_DA*)da->data;
735: if (lx) *lx = dd->lx;
736: if (ly) *ly = dd->ly;
737: if (lz) *lz = dd->lz;
738: return(0);
739: }
741: /*@
742: DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
744: Logically Collective on DMDA
746: Input Parameters:
747: + da - the DMDA object
748: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
749: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
750: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
752: Options Database:
753: + -da_refine_x - refinement ratio in x direction
754: . -da_refine_y - refinement ratio in y direction
755: - -da_refine_z - refinement ratio in z direction
757: Level: intermediate
759: Notes: Pass PETSC_IGNORE to leave a value unchanged
761: .seealso: DMRefine(), DMDAGetRefinementFactor()
762: @*/
763: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
764: {
765: DM_DA *dd = (DM_DA*)da->data;
773: if (refine_x > 0) dd->refine_x = refine_x;
774: if (refine_y > 0) dd->refine_y = refine_y;
775: if (refine_z > 0) dd->refine_z = refine_z;
776: return(0);
777: }
779: /*@C
780: DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
782: Not Collective
784: Input Parameter:
785: . da - the DMDA object
787: Output Parameters:
788: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
789: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
790: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
792: Level: intermediate
794: Notes: Pass NULL for values you do not need
796: .seealso: DMRefine(), DMDASetRefinementFactor()
797: @*/
798: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
799: {
800: DM_DA *dd = (DM_DA*)da->data;
804: if (refine_x) *refine_x = dd->refine_x;
805: if (refine_y) *refine_y = dd->refine_y;
806: if (refine_z) *refine_z = dd->refine_z;
807: return(0);
808: }
810: /*@C
811: DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
813: Logically Collective on DMDA
815: Input Parameters:
816: + da - the DMDA object
817: - f - the function that allocates the matrix for that specific DMDA
819: Level: developer
821: Notes: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
822: the diagonal and off-diagonal blocks of the matrix
824: Not supported from Fortran
826: .seealso: DMCreateMatrix(), DMDASetBlockFills()
827: @*/
828: PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
829: {
832: da->ops->creatematrix = f;
833: return(0);
834: }
836: /*
837: Creates "balanced" ownership ranges after refinement, constrained by the need for the
838: fine grid boundaries to fall within one stencil width of the coarse partition.
840: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
841: */
842: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
843: {
844: PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
848: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
849: if (ratio == 1) {
850: PetscMemcpy(lf,lc,m*sizeof(lc[0]));
851: return(0);
852: }
853: for (i=0; i<m; i++) totalc += lc[i];
854: remaining = (!periodic) + ratio * (totalc - (!periodic));
855: for (i=0; i<m; i++) {
856: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
857: if (i == m-1) lf[i] = want;
858: else {
859: const PetscInt nextc = startc + lc[i];
860: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
861: * coarse stencil width of the first coarse node in the next subdomain. */
862: while ((startf+want)/ratio < nextc - stencil_width) want++;
863: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
864: * coarse stencil width of the last coarse node in the current subdomain. */
865: while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
866: /* Make sure all constraints are satisfied */
867: if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
868: || ((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");
869: }
870: lf[i] = want;
871: startc += lc[i];
872: startf += lf[i];
873: remaining -= lf[i];
874: }
875: return(0);
876: }
878: /*
879: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
880: fine grid boundaries to fall within one stencil width of the coarse partition.
882: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
883: */
884: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
885: {
886: PetscInt i,totalf,remaining,startc,startf;
890: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
891: if (ratio == 1) {
892: PetscMemcpy(lc,lf,m*sizeof(lf[0]));
893: return(0);
894: }
895: for (i=0,totalf=0; i<m; i++) totalf += lf[i];
896: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
897: for (i=0,startc=0,startf=0; i<m; i++) {
898: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
899: if (i == m-1) lc[i] = want;
900: else {
901: const PetscInt nextf = startf+lf[i];
902: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
903: * node is within one stencil width. */
904: while (nextf/ratio < startc+want-stencil_width) want--;
905: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
906: * fine node is within one stencil width. */
907: while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
908: if (want < 0 || want > remaining
909: || (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");
910: }
911: lc[i] = want;
912: startc += lc[i];
913: startf += lf[i];
914: remaining -= lc[i];
915: }
916: return(0);
917: }
919: PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
920: {
922: PetscInt M,N,P,i,dim;
923: DM da2;
924: DM_DA *dd = (DM_DA*)da->data,*dd2;
930: DMGetDimension(da, &dim);
931: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
932: M = dd->refine_x*dd->M;
933: } else {
934: M = 1 + dd->refine_x*(dd->M - 1);
935: }
936: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
937: if (dim > 1) {
938: N = dd->refine_y*dd->N;
939: } else {
940: N = 1;
941: }
942: } else {
943: N = 1 + dd->refine_y*(dd->N - 1);
944: }
945: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
946: if (dim > 2) {
947: P = dd->refine_z*dd->P;
948: } else {
949: P = 1;
950: }
951: } else {
952: P = 1 + dd->refine_z*(dd->P - 1);
953: }
954: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
955: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
956: DMSetDimension(da2,dim);
957: DMDASetSizes(da2,M,N,P);
958: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
959: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
960: DMDASetDof(da2,dd->w);
961: DMDASetStencilType(da2,dd->stencil_type);
962: DMDASetStencilWidth(da2,dd->s);
963: if (dim == 3) {
964: PetscInt *lx,*ly,*lz;
965: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
966: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
967: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
968: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
969: DMDASetOwnershipRanges(da2,lx,ly,lz);
970: PetscFree3(lx,ly,lz);
971: } else if (dim == 2) {
972: PetscInt *lx,*ly;
973: PetscMalloc2(dd->m,&lx,dd->n,&ly);
974: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
975: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
976: DMDASetOwnershipRanges(da2,lx,ly,NULL);
977: PetscFree2(lx,ly);
978: } else if (dim == 1) {
979: PetscInt *lx;
980: PetscMalloc1(dd->m,&lx);
981: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
982: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
983: PetscFree(lx);
984: }
985: dd2 = (DM_DA*)da2->data;
987: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
988: da2->ops->creatematrix = da->ops->creatematrix;
989: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
990: da2->ops->getcoloring = da->ops->getcoloring;
991: dd2->interptype = dd->interptype;
993: /* copy fill information if given */
994: if (dd->dfill) {
995: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
996: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
997: }
998: if (dd->ofill) {
999: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1000: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1001: }
1002: /* copy the refine information */
1003: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
1004: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
1005: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
1007: if (dd->refine_z_hier) {
1008: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) {
1009: dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
1010: }
1011: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) {
1012: dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
1013: }
1014: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1015: PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1016: PetscMemcpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n*sizeof(PetscInt));
1017: }
1018: if (dd->refine_y_hier) {
1019: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) {
1020: dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
1021: }
1022: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) {
1023: dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
1024: }
1025: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1026: PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1027: PetscMemcpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n*sizeof(PetscInt));
1028: }
1029: if (dd->refine_x_hier) {
1030: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) {
1031: dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
1032: }
1033: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) {
1034: dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
1035: }
1036: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1037: PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1038: PetscMemcpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n*sizeof(PetscInt));
1039: }
1040:
1042: /* copy vector type information */
1043: DMSetVecType(da2,da->vectype);
1045: dd2->lf = dd->lf;
1046: dd2->lj = dd->lj;
1048: da2->leveldown = da->leveldown;
1049: da2->levelup = da->levelup + 1;
1051: DMSetUp(da2);
1053: /* interpolate coordinates if they are set on the coarse grid */
1054: if (da->coordinates) {
1055: DM cdaf,cdac;
1056: Vec coordsc,coordsf;
1057: Mat II;
1059: DMGetCoordinateDM(da,&cdac);
1060: DMGetCoordinates(da,&coordsc);
1061: DMGetCoordinateDM(da2,&cdaf);
1062: /* force creation of the coordinate vector */
1063: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1064: DMGetCoordinates(da2,&coordsf);
1065: DMCreateInterpolation(cdac,cdaf,&II,NULL);
1066: MatInterpolate(II,coordsc,coordsf);
1067: MatDestroy(&II);
1068: }
1070: for (i=0; i<da->bs; i++) {
1071: const char *fieldname;
1072: DMDAGetFieldName(da,i,&fieldname);
1073: DMDASetFieldName(da2,i,fieldname);
1074: }
1076: *daref = da2;
1077: return(0);
1078: }
1081: PetscErrorCode DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1082: {
1084: PetscInt M,N,P,i,dim;
1085: DM da2;
1086: DM_DA *dd = (DM_DA*)da->data,*dd2;
1092: DMGetDimension(da, &dim);
1093: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1094: M = dd->M / dd->coarsen_x;
1095: } else {
1096: M = 1 + (dd->M - 1) / dd->coarsen_x;
1097: }
1098: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1099: if (dim > 1) {
1100: N = dd->N / dd->coarsen_y;
1101: } else {
1102: N = 1;
1103: }
1104: } else {
1105: N = 1 + (dd->N - 1) / dd->coarsen_y;
1106: }
1107: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1108: if (dim > 2) {
1109: P = dd->P / dd->coarsen_z;
1110: } else {
1111: P = 1;
1112: }
1113: } else {
1114: P = 1 + (dd->P - 1) / dd->coarsen_z;
1115: }
1116: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1117: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1118: DMSetDimension(da2,dim);
1119: DMDASetSizes(da2,M,N,P);
1120: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1121: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1122: DMDASetDof(da2,dd->w);
1123: DMDASetStencilType(da2,dd->stencil_type);
1124: DMDASetStencilWidth(da2,dd->s);
1125: if (dim == 3) {
1126: PetscInt *lx,*ly,*lz;
1127: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1128: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1129: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1130: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1131: DMDASetOwnershipRanges(da2,lx,ly,lz);
1132: PetscFree3(lx,ly,lz);
1133: } else if (dim == 2) {
1134: PetscInt *lx,*ly;
1135: PetscMalloc2(dd->m,&lx,dd->n,&ly);
1136: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1137: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1138: DMDASetOwnershipRanges(da2,lx,ly,NULL);
1139: PetscFree2(lx,ly);
1140: } else if (dim == 1) {
1141: PetscInt *lx;
1142: PetscMalloc1(dd->m,&lx);
1143: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1144: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1145: PetscFree(lx);
1146: }
1147: dd2 = (DM_DA*)da2->data;
1149: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1150: /* da2->ops->createinterpolation = da->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1151: da2->ops->creatematrix = da->ops->creatematrix;
1152: da2->ops->getcoloring = da->ops->getcoloring;
1153: dd2->interptype = dd->interptype;
1155: /* copy fill information if given */
1156: if (dd->dfill) {
1157: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1158: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1159: }
1160: if (dd->ofill) {
1161: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1162: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1163: }
1164: /* copy the refine information */
1165: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1166: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1167: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1169: if (dd->refine_z_hier) {
1170: if (da->levelup - da->leveldown -1 > -1 && da->levelup - da->leveldown - 1< dd->refine_z_hier_n) {
1171: dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown - 1];
1172: }
1173: if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_z_hier_n) {
1174: dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown - 2];
1175: }
1176: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1177: PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1178: PetscMemcpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n*sizeof(PetscInt));
1179: }
1180: if (dd->refine_y_hier) {
1181: if (da->levelup - da->leveldown - 1 > -1 && da->levelup - da->leveldown - 1< dd->refine_y_hier_n) {
1182: dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown - 1];
1183: }
1184: if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_y_hier_n) {
1185: dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown - 2];
1186: }
1187: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1188: PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1189: PetscMemcpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n*sizeof(PetscInt));
1190: }
1191: if (dd->refine_x_hier) {
1192: if (da->levelup - da->leveldown - 1 > -1 && da->levelup - da->leveldown - 1 < dd->refine_x_hier_n) {
1193: dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown - 1];
1194: }
1195: if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_x_hier_n) {
1196: dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown - 2];
1197: }
1198: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1199: PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1200: PetscMemcpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n*sizeof(PetscInt));
1201: }
1203: /* copy vector type information */
1204: DMSetVecType(da2,da->vectype);
1206: dd2->lf = dd->lf;
1207: dd2->lj = dd->lj;
1209: da2->leveldown = da->leveldown + 1;
1210: da2->levelup = da->levelup;
1212: DMSetUp(da2);
1214: /* inject coordinates if they are set on the fine grid */
1215: if (da->coordinates) {
1216: DM cdaf,cdac;
1217: Vec coordsc,coordsf;
1218: Mat inject;
1219: VecScatter vscat;
1221: DMGetCoordinateDM(da,&cdaf);
1222: DMGetCoordinates(da,&coordsf);
1223: DMGetCoordinateDM(da2,&cdac);
1224: /* force creation of the coordinate vector */
1225: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1226: DMGetCoordinates(da2,&coordsc);
1228: DMCreateInjection(cdac,cdaf,&inject);
1229: MatScatterGetVecScatter(inject,&vscat);
1230: VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1231: VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1232: MatDestroy(&inject);
1233: }
1235: for (i=0; i<da->bs; i++) {
1236: const char *fieldname;
1237: DMDAGetFieldName(da,i,&fieldname);
1238: DMDASetFieldName(da2,i,fieldname);
1239: }
1241: *daref = da2;
1242: return(0);
1243: }
1245: PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1246: {
1248: PetscInt i,n,*refx,*refy,*refz;
1252: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1253: if (nlevels == 0) return(0);
1256: /* Get refinement factors, defaults taken from the coarse DMDA */
1257: PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1258: for (i=0; i<nlevels; i++) {
1259: DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1260: }
1261: n = nlevels;
1262: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1263: n = nlevels;
1264: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1265: n = nlevels;
1266: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);
1268: DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1269: DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1270: for (i=1; i<nlevels; i++) {
1271: DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1272: DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1273: }
1274: PetscFree3(refx,refy,refz);
1275: return(0);
1276: }
1278: PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1279: {
1281: PetscInt i;
1285: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1286: if (nlevels == 0) return(0);
1288: DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1289: for (i=1; i<nlevels; i++) {
1290: DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1291: }
1292: return(0);
1293: }
1295: #include <petscgll.h>
1297: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm,PetscGLL *gll)
1298: {
1300: PetscInt i,j,n = gll->n,xs,xn,q;
1301: PetscScalar *xx;
1302: PetscReal h;
1303: Vec x;
1304: DM_DA *da = (DM_DA*)dm->data;
1307: if (da->bx != DM_BOUNDARY_PERIODIC) {
1308: DMDAGetInfo(dm,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1309: q = (q-1)/(n-1); /* number of spectral elements */
1310: h = 2.0/q;
1311: DMDAGetCorners(dm,&xs,NULL,NULL,&xn,NULL,NULL);
1312: xs = xs/(n-1);
1313: xn = xn/(n-1);
1314: DMDASetUniformCoordinates(dm,-1.,1.,0.,0.,0.,0.);
1315: DMGetCoordinates(dm,&x);
1316: DMDAVecGetArray(dm,x,&xx);
1318: /* loop over local spectral elements */
1319: for (j=xs; j<xs+xn; j++) {
1320: /*
1321: Except for the first process, each process starts on the second GLL point of the first element on that process
1322: */
1323: for (i= (j == xs && xs > 0)? 1 : 0; i<n; i++) {
1324: xx[j*(n-1) + i] = -1.0 + h*j + h*(gll->nodes[i]+1.0)/2.;
1325: }
1326: }
1327: DMDAVecRestoreArray(dm,x,&xx);
1328: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for periodic");
1329: return(0);
1330: }
1332: /*@
1334: 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
1336: Collective on DM
1338: Input Parameters:
1339: + da - the DMDA object
1340: - gll - the GLL object
1342: Notes: the parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1343: on each process much be divisible by the number of GLL elements needed per process. This depends on whether the DM is
1344: periodic or not.
1346: Level: advanced
1348: .seealso: DMDACreate(), PetscGLLCreate(), DMGetCoordinates()
1349: @*/
1350: PetscErrorCode DMDASetGLLCoordinates(DM da,PetscGLL *gll)
1351: {
1355: if (da->dim == 1) {
1356: DMDASetGLLCoordinates_1d(da,gll);
1357: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for 2 or 3d");
1358: return(0);
1359: }