Actual source code: da.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/dmdaimpl.h>
3: /*@
4: DMDASetSizes - Sets the number of grid points in the three dimensional directions
6: Logically Collective on da
8: Input Parameters:
9: + da - the DMDA
10: . M - the global X size
11: . N - the global Y size
12: - P - the global Z size
14: Level: intermediate
16: Developer Notes:
17: Since the dimension may not yet have been set the code cannot error check for non-positive Y and Z number of grid points
19: .seealso: PetscSplitOwnership()
20: @*/
21: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
22: {
23: DM_DA *dd = (DM_DA*)da->data;
30: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
31: if (M < 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in X direction must be positive");
32: if (N < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in Y direction must be positive");
33: if (P < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in Z direction must be positive");
35: dd->M = M;
36: dd->N = N;
37: dd->P = P;
38: return(0);
39: }
41: /*@
42: DMDASetNumProcs - Sets the number of processes in each dimension
44: Logically Collective on da
46: Input Parameters:
47: + da - the DMDA
48: . m - the number of X procs (or PETSC_DECIDE)
49: . n - the number of Y procs (or PETSC_DECIDE)
50: - p - the number of Z procs (or PETSC_DECIDE)
52: Level: intermediate
54: .seealso: DMDASetSizes(), PetscSplitOwnership()
55: @*/
56: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
57: {
58: DM_DA *dd = (DM_DA*)da->data;
66: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
67: dd->m = m;
68: dd->n = n;
69: dd->p = p;
70: if (da->dim == 2) {
71: PetscMPIInt size;
72: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
73: if ((dd->m > 0) && (dd->n < 0)) {
74: dd->n = size/dd->m;
75: 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);
76: }
77: if ((dd->n > 0) && (dd->m < 0)) {
78: dd->m = size/dd->n;
79: 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);
80: }
81: }
82: return(0);
83: }
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: .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: .seealso: DMDAGetDof(), DMDACreate(), DMDestroy(), DMDA
126: @*/
127: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
128: {
129: DM_DA *dd = (DM_DA*)da->data;
134: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
135: dd->w = dof;
136: da->bs = dof;
137: return(0);
138: }
140: /*@
141: DMDAGetDof - Gets the number of degrees of freedom per vertex
143: Not collective
145: Input Parameter:
146: . da - The DMDA
148: Output Parameter:
149: . dof - Number of degrees of freedom
151: Level: intermediate
153: .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
154: @*/
155: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
156: {
157: DM_DA *dd = (DM_DA *) da->data;
162: *dof = dd->w;
163: return(0);
164: }
166: /*@
167: DMDAGetOverlap - Gets the size of the per-processor overlap.
169: Not collective
171: Input Parameters:
172: . da - The DMDA
174: Output Parameters:
175: + x - Overlap in the x direction
176: . y - Overlap in the y direction
177: - z - Overlap in the z direction
179: Level: intermediate
181: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
182: @*/
183: PetscErrorCode DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
184: {
185: DM_DA *dd = (DM_DA*)da->data;
189: if (x) *x = dd->xol;
190: if (y) *y = dd->yol;
191: if (z) *z = dd->zol;
192: return(0);
193: }
195: /*@
196: DMDASetOverlap - Sets the size of the per-processor overlap.
198: Not collective
200: Input Parameters:
201: + da - The DMDA
202: . x - Overlap in the x direction
203: . y - Overlap in the y direction
204: - z - Overlap in the z direction
206: Level: intermediate
208: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
209: @*/
210: PetscErrorCode DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
211: {
212: DM_DA *dd = (DM_DA*)da->data;
219: dd->xol = x;
220: dd->yol = y;
221: dd->zol = z;
222: return(0);
223: }
226: /*@
227: DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.
229: Not collective
231: Input Parameters:
232: . da - The DMDA
234: Output Parameters:
235: . Nsub - Number of local subdomains created upon decomposition
237: Level: intermediate
239: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
240: @*/
241: PetscErrorCode DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
242: {
243: DM_DA *dd = (DM_DA*)da->data;
247: if (Nsub) *Nsub = dd->Nsub;
248: return(0);
249: }
251: /*@
252: DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.
254: Not collective
256: Input Parameters:
257: + da - The DMDA
258: - Nsub - The number of local subdomains requested
260: Level: intermediate
262: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
263: @*/
264: PetscErrorCode DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
265: {
266: DM_DA *dd = (DM_DA*)da->data;
271: dd->Nsub = Nsub;
272: return(0);
273: }
275: /*@
276: DMDASetOffset - Sets the index offset of the DA.
278: Collective on DA
280: Input Parameter:
281: + da - The DMDA
282: . xo - The offset in the x direction
283: . yo - The offset in the y direction
284: - zo - The offset in the z direction
286: Level: intermediate
288: Notes:
289: This is used primarily to overlap a computation on a local DA with that on a global DA without
290: changing boundary conditions or subdomain features that depend upon the global offsets.
292: .seealso: DMDAGetOffset(), DMDAVecGetArray()
293: @*/
294: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
295: {
297: DM_DA *dd = (DM_DA*)da->data;
307: dd->xo = xo;
308: dd->yo = yo;
309: dd->zo = zo;
310: dd->Mo = Mo;
311: dd->No = No;
312: dd->Po = Po;
314: if (da->coordinateDM) {
315: DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
316: }
317: return(0);
318: }
320: /*@
321: DMDAGetOffset - Gets the index offset of the DA.
323: Not collective
325: Input Parameter:
326: . da - The DMDA
328: Output Parameters:
329: + xo - The offset in the x direction
330: . yo - The offset in the y direction
331: . zo - The offset in the z direction
332: . Mo - The global size in the x direction
333: . No - The global size in the y direction
334: - Po - The global size in the z direction
336: Level: intermediate
338: .seealso: DMDASetOffset(), DMDAVecGetArray()
339: @*/
340: PetscErrorCode DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
341: {
342: DM_DA *dd = (DM_DA*)da->data;
346: if (xo) *xo = dd->xo;
347: if (yo) *yo = dd->yo;
348: if (zo) *zo = dd->zo;
349: if (Mo) *Mo = dd->Mo;
350: if (No) *No = dd->No;
351: if (Po) *Po = dd->Po;
352: return(0);
353: }
355: /*@
356: DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.
358: Not collective
360: Input Parameter:
361: . da - The DMDA
363: Output Parameters:
364: + xs - The start of the region in x
365: . ys - The start of the region in y
366: . zs - The start of the region in z
367: . xs - The size of the region in x
368: . ys - The size of the region in y
369: - zs - The size of the region in z
371: Level: intermediate
373: .seealso: DMDAGetOffset(), DMDAVecGetArray()
374: @*/
375: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
376: {
377: DM_DA *dd = (DM_DA*)da->data;
381: if (xs) *xs = dd->nonxs;
382: if (ys) *ys = dd->nonys;
383: if (zs) *zs = dd->nonzs;
384: if (xm) *xm = dd->nonxm;
385: if (ym) *ym = dd->nonym;
386: if (zm) *zm = dd->nonzm;
387: return(0);
388: }
391: /*@
392: DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.
394: Collective on DA
396: Input Parameter:
397: + da - The DMDA
398: . xs - The start of the region in x
399: . ys - The start of the region in y
400: . zs - The start of the region in z
401: . xs - The size of the region in x
402: . ys - The size of the region in y
403: - zs - The size of the region in z
405: Level: intermediate
407: .seealso: DMDAGetOffset(), DMDAVecGetArray()
408: @*/
409: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
410: {
411: DM_DA *dd = (DM_DA*)da->data;
421: dd->nonxs = xs;
422: dd->nonys = ys;
423: dd->nonzs = zs;
424: dd->nonxm = xm;
425: dd->nonym = ym;
426: dd->nonzm = zm;
428: return(0);
429: }
431: /*@
432: DMDASetStencilType - Sets the type of the communication stencil
434: Logically Collective on da
436: Input Parameter:
437: + da - The DMDA
438: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
440: Level: intermediate
442: .seealso: DMDACreate(), DMDestroy(), DMDA
443: @*/
444: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
445: {
446: DM_DA *dd = (DM_DA*)da->data;
451: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
452: dd->stencil_type = stype;
453: return(0);
454: }
456: /*@
457: DMDAGetStencilType - Gets the type of the communication stencil
459: Not collective
461: Input Parameter:
462: . da - The DMDA
464: Output Parameter:
465: . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
467: Level: intermediate
469: .seealso: DMDACreate(), DMDestroy(), DMDA
470: @*/
471: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
472: {
473: DM_DA *dd = (DM_DA*)da->data;
478: *stype = dd->stencil_type;
479: return(0);
480: }
482: /*@
483: DMDASetStencilWidth - Sets the width of the communication stencil
485: Logically Collective on da
487: Input Parameter:
488: + da - The DMDA
489: - width - The stencil width
491: Level: intermediate
493: .seealso: DMDACreate(), DMDestroy(), DMDA
494: @*/
495: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
496: {
497: DM_DA *dd = (DM_DA*)da->data;
502: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
503: dd->s = width;
504: return(0);
505: }
507: /*@
508: DMDAGetStencilWidth - Gets the width of the communication stencil
510: Not collective
512: Input Parameter:
513: . da - The DMDA
515: Output Parameter:
516: . width - The stencil width
518: Level: intermediate
520: .seealso: DMDACreate(), DMDestroy(), DMDA
521: @*/
522: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
523: {
524: DM_DA *dd = (DM_DA *) da->data;
529: *width = dd->s;
530: return(0);
531: }
533: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
534: {
535: PetscInt i,sum;
538: if (M < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
539: for (i=sum=0; i<m; i++) sum += lx[i];
540: if (sum != M) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
541: return(0);
542: }
544: /*@
545: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
547: Logically Collective on da
549: Input Parameter:
550: + da - The DMDA
551: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
552: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
553: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.
555: Level: intermediate
557: Note: these numbers are NOT multiplied by the number of dof per node.
559: .seealso: DMDACreate(), DMDestroy(), DMDA
560: @*/
561: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
562: {
564: DM_DA *dd = (DM_DA*)da->data;
568: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
569: if (lx) {
570: if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
571: DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
572: if (!dd->lx) {
573: PetscMalloc1(dd->m, &dd->lx);
574: }
575: PetscArraycpy(dd->lx, lx, dd->m);
576: }
577: if (ly) {
578: if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
579: DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
580: if (!dd->ly) {
581: PetscMalloc1(dd->n, &dd->ly);
582: }
583: PetscArraycpy(dd->ly, ly, dd->n);
584: }
585: if (lz) {
586: if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
587: DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
588: if (!dd->lz) {
589: PetscMalloc1(dd->p, &dd->lz);
590: }
591: PetscArraycpy(dd->lz, lz, dd->p);
592: }
593: return(0);
594: }
596: /*@
597: DMDASetInterpolationType - Sets the type of interpolation that will be
598: returned by DMCreateInterpolation()
600: Logically Collective on da
602: Input Parameter:
603: + da - initial distributed array
604: - ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
606: Level: intermediate
608: Notes:
609: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
611: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
612: @*/
613: PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
614: {
615: DM_DA *dd = (DM_DA*)da->data;
620: dd->interptype = ctype;
621: return(0);
622: }
624: /*@
625: DMDAGetInterpolationType - Gets the type of interpolation that will be
626: used by DMCreateInterpolation()
628: Not Collective
630: Input Parameter:
631: . da - distributed array
633: Output Parameter:
634: . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
636: Level: intermediate
638: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
639: @*/
640: PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
641: {
642: DM_DA *dd = (DM_DA*)da->data;
647: *ctype = dd->interptype;
648: return(0);
649: }
651: /*@C
652: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
653: processes neighbors.
655: Not Collective
657: Input Parameter:
658: . da - the DMDA object
660: Output Parameters:
661: . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
662: this process itself is in the list
664: Notes:
665: In 2d the array is of length 9, in 3d of length 27
666: Not supported in 1d
667: Do not free the array, it is freed when the DMDA is destroyed.
669: Fortran Notes:
670: In fortran you must pass in an array of the appropriate length.
672: Level: intermediate
674: @*/
675: PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
676: {
677: DM_DA *dd = (DM_DA*)da->data;
681: *ranks = dd->neighbors;
682: return(0);
683: }
685: /*@C
686: DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
688: Not Collective
690: Input Parameter:
691: . da - the DMDA object
693: Output Parameter:
694: + lx - ownership along x direction (optional)
695: . ly - ownership along y direction (optional)
696: - lz - ownership along z direction (optional)
698: Level: intermediate
700: Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
702: In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
703: eighth arguments from DMDAGetInfo()
705: In C you should not free these arrays, nor change the values in them. They will only have valid values while the
706: DMDA they came from still exists (has not been destroyed).
708: These numbers are NOT multiplied by the number of dof per node.
710: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
711: @*/
712: PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
713: {
714: DM_DA *dd = (DM_DA*)da->data;
718: if (lx) *lx = dd->lx;
719: if (ly) *ly = dd->ly;
720: if (lz) *lz = dd->lz;
721: return(0);
722: }
724: /*@
725: DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
727: Logically Collective on da
729: Input Parameters:
730: + da - the DMDA object
731: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
732: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
733: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
735: Options Database:
736: + -da_refine_x refine_x - refinement ratio in x direction
737: . -da_refine_y rafine_y - refinement ratio in y direction
738: . -da_refine_z refine_z - refinement ratio in z direction
739: - -da_refine <n> - refine the DMDA object n times when it is created.
741: Level: intermediate
743: Notes:
744: Pass PETSC_IGNORE to leave a value unchanged
746: .seealso: DMRefine(), DMDAGetRefinementFactor()
747: @*/
748: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
749: {
750: DM_DA *dd = (DM_DA*)da->data;
758: if (refine_x > 0) dd->refine_x = refine_x;
759: if (refine_y > 0) dd->refine_y = refine_y;
760: if (refine_z > 0) dd->refine_z = refine_z;
761: return(0);
762: }
764: /*@C
765: DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
767: Not Collective
769: Input Parameter:
770: . da - the DMDA object
772: Output Parameters:
773: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
774: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
775: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
777: Level: intermediate
779: Notes:
780: Pass NULL for values you do not need
782: .seealso: DMRefine(), DMDASetRefinementFactor()
783: @*/
784: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
785: {
786: DM_DA *dd = (DM_DA*)da->data;
790: if (refine_x) *refine_x = dd->refine_x;
791: if (refine_y) *refine_y = dd->refine_y;
792: if (refine_z) *refine_z = dd->refine_z;
793: return(0);
794: }
796: /*@C
797: DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
799: Logically Collective on da
801: Input Parameters:
802: + da - the DMDA object
803: - f - the function that allocates the matrix for that specific DMDA
805: Level: developer
807: Notes:
808: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
809: the diagonal and off-diagonal blocks of the matrix
811: Not supported from Fortran
813: .seealso: DMCreateMatrix(), DMDASetBlockFills()
814: @*/
815: PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
816: {
819: da->ops->creatematrix = f;
820: return(0);
821: }
823: /*
824: Creates "balanced" ownership ranges after refinement, constrained by the need for the
825: fine grid boundaries to fall within one stencil width of the coarse partition.
827: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
828: */
829: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
830: {
831: PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
835: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
836: if (ratio == 1) {
837: PetscArraycpy(lf,lc,m);
838: return(0);
839: }
840: for (i=0; i<m; i++) totalc += lc[i];
841: remaining = (!periodic) + ratio * (totalc - (!periodic));
842: for (i=0; i<m; i++) {
843: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
844: if (i == m-1) lf[i] = want;
845: else {
846: const PetscInt nextc = startc + lc[i];
847: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
848: * coarse stencil width of the first coarse node in the next subdomain. */
849: while ((startf+want)/ratio < nextc - stencil_width) want++;
850: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
851: * coarse stencil width of the last coarse node in the current subdomain. */
852: while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
853: /* Make sure all constraints are satisfied */
854: if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
855: || ((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");
856: }
857: lf[i] = want;
858: startc += lc[i];
859: startf += lf[i];
860: remaining -= lf[i];
861: }
862: return(0);
863: }
865: /*
866: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
867: fine grid boundaries to fall within one stencil width of the coarse partition.
869: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
870: */
871: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
872: {
873: PetscInt i,totalf,remaining,startc,startf;
877: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
878: if (ratio == 1) {
879: PetscArraycpy(lc,lf,m);
880: return(0);
881: }
882: for (i=0,totalf=0; i<m; i++) totalf += lf[i];
883: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
884: for (i=0,startc=0,startf=0; i<m; i++) {
885: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
886: if (i == m-1) lc[i] = want;
887: else {
888: const PetscInt nextf = startf+lf[i];
889: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
890: * node is within one stencil width. */
891: while (nextf/ratio < startc+want-stencil_width) want--;
892: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
893: * fine node is within one stencil width. */
894: while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
895: if (want < 0 || want > remaining
896: || (nextf/ratio < startc+want-stencil_width) || ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible coarsened ownership range");
897: }
898: lc[i] = want;
899: startc += lc[i];
900: startf += lf[i];
901: remaining -= lc[i];
902: }
903: return(0);
904: }
906: PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
907: {
909: PetscInt M,N,P,i,dim;
910: DM da2;
911: DM_DA *dd = (DM_DA*)da->data,*dd2;
917: DMGetDimension(da, &dim);
918: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
919: M = dd->refine_x*dd->M;
920: } else {
921: M = 1 + dd->refine_x*(dd->M - 1);
922: }
923: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
924: if (dim > 1) {
925: N = dd->refine_y*dd->N;
926: } else {
927: N = 1;
928: }
929: } else {
930: N = 1 + dd->refine_y*(dd->N - 1);
931: }
932: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
933: if (dim > 2) {
934: P = dd->refine_z*dd->P;
935: } else {
936: P = 1;
937: }
938: } else {
939: P = 1 + dd->refine_z*(dd->P - 1);
940: }
941: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
942: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
943: DMSetDimension(da2,dim);
944: DMDASetSizes(da2,M,N,P);
945: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
946: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
947: DMDASetDof(da2,dd->w);
948: DMDASetStencilType(da2,dd->stencil_type);
949: DMDASetStencilWidth(da2,dd->s);
950: if (dim == 3) {
951: PetscInt *lx,*ly,*lz;
952: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
953: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
954: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
955: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
956: DMDASetOwnershipRanges(da2,lx,ly,lz);
957: PetscFree3(lx,ly,lz);
958: } else if (dim == 2) {
959: PetscInt *lx,*ly;
960: PetscMalloc2(dd->m,&lx,dd->n,&ly);
961: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
962: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
963: DMDASetOwnershipRanges(da2,lx,ly,NULL);
964: PetscFree2(lx,ly);
965: } else if (dim == 1) {
966: PetscInt *lx;
967: PetscMalloc1(dd->m,&lx);
968: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
969: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
970: PetscFree(lx);
971: }
972: dd2 = (DM_DA*)da2->data;
974: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
975: da2->ops->creatematrix = da->ops->creatematrix;
976: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
977: da2->ops->getcoloring = da->ops->getcoloring;
978: dd2->interptype = dd->interptype;
980: /* copy fill information if given */
981: if (dd->dfill) {
982: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
983: PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);
984: }
985: if (dd->ofill) {
986: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
987: PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);
988: }
989: /* copy the refine information */
990: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
991: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
992: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
994: if (dd->refine_z_hier) {
995: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) {
996: dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
997: }
998: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) {
999: dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
1000: }
1001: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1002: PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1003: PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);
1004: }
1005: if (dd->refine_y_hier) {
1006: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) {
1007: dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
1008: }
1009: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) {
1010: dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
1011: }
1012: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1013: PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1014: PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);
1015: }
1016: if (dd->refine_x_hier) {
1017: if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) {
1018: dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
1019: }
1020: if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) {
1021: dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
1022: }
1023: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1024: PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1025: PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);
1026: }
1029: /* copy vector type information */
1030: DMSetVecType(da2,da->vectype);
1032: dd2->lf = dd->lf;
1033: dd2->lj = dd->lj;
1035: da2->leveldown = da->leveldown;
1036: da2->levelup = da->levelup + 1;
1038: DMSetUp(da2);
1040: /* interpolate coordinates if they are set on the coarse grid */
1041: if (da->coordinates) {
1042: DM cdaf,cdac;
1043: Vec coordsc,coordsf;
1044: Mat II;
1046: DMGetCoordinateDM(da,&cdac);
1047: DMGetCoordinates(da,&coordsc);
1048: DMGetCoordinateDM(da2,&cdaf);
1049: /* force creation of the coordinate vector */
1050: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1051: DMGetCoordinates(da2,&coordsf);
1052: DMCreateInterpolation(cdac,cdaf,&II,NULL);
1053: MatInterpolate(II,coordsc,coordsf);
1054: MatDestroy(&II);
1055: }
1057: for (i=0; i<da->bs; i++) {
1058: const char *fieldname;
1059: DMDAGetFieldName(da,i,&fieldname);
1060: DMDASetFieldName(da2,i,fieldname);
1061: }
1063: *daref = da2;
1064: return(0);
1065: }
1068: PetscErrorCode DMCoarsen_DA(DM dmf, MPI_Comm comm,DM *dmc)
1069: {
1071: PetscInt M,N,P,i,dim;
1072: DM dmc2;
1073: DM_DA *dd = (DM_DA*)dmf->data,*dd2;
1079: DMGetDimension(dmf, &dim);
1080: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1081: M = dd->M / dd->coarsen_x;
1082: } else {
1083: M = 1 + (dd->M - 1) / dd->coarsen_x;
1084: }
1085: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1086: if (dim > 1) {
1087: N = dd->N / dd->coarsen_y;
1088: } else {
1089: N = 1;
1090: }
1091: } else {
1092: N = 1 + (dd->N - 1) / dd->coarsen_y;
1093: }
1094: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1095: if (dim > 2) {
1096: P = dd->P / dd->coarsen_z;
1097: } else {
1098: P = 1;
1099: }
1100: } else {
1101: P = 1 + (dd->P - 1) / dd->coarsen_z;
1102: }
1103: DMDACreate(PetscObjectComm((PetscObject)dmf),&dmc2);
1104: DMSetOptionsPrefix(dmc2,((PetscObject)dmf)->prefix);
1105: DMSetDimension(dmc2,dim);
1106: DMDASetSizes(dmc2,M,N,P);
1107: DMDASetNumProcs(dmc2,dd->m,dd->n,dd->p);
1108: DMDASetBoundaryType(dmc2,dd->bx,dd->by,dd->bz);
1109: DMDASetDof(dmc2,dd->w);
1110: DMDASetStencilType(dmc2,dd->stencil_type);
1111: DMDASetStencilWidth(dmc2,dd->s);
1112: if (dim == 3) {
1113: PetscInt *lx,*ly,*lz;
1114: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1115: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1116: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1117: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1118: DMDASetOwnershipRanges(dmc2,lx,ly,lz);
1119: PetscFree3(lx,ly,lz);
1120: } else if (dim == 2) {
1121: PetscInt *lx,*ly;
1122: PetscMalloc2(dd->m,&lx,dd->n,&ly);
1123: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1124: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1125: DMDASetOwnershipRanges(dmc2,lx,ly,NULL);
1126: PetscFree2(lx,ly);
1127: } else if (dim == 1) {
1128: PetscInt *lx;
1129: PetscMalloc1(dd->m,&lx);
1130: DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1131: DMDASetOwnershipRanges(dmc2,lx,NULL,NULL);
1132: PetscFree(lx);
1133: }
1134: dd2 = (DM_DA*)dmc2->data;
1136: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1137: /* dmc2->ops->createinterpolation = dmf->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1138: dmc2->ops->creatematrix = dmf->ops->creatematrix;
1139: dmc2->ops->getcoloring = dmf->ops->getcoloring;
1140: dd2->interptype = dd->interptype;
1142: /* copy fill information if given */
1143: if (dd->dfill) {
1144: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1145: PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);
1146: }
1147: if (dd->ofill) {
1148: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1149: PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);
1150: }
1151: /* copy the refine information */
1152: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1153: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1154: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1156: if (dd->refine_z_hier) {
1157: if (dmf->levelup - dmf->leveldown -1 > -1 && dmf->levelup - dmf->leveldown - 1< dd->refine_z_hier_n) {
1158: dd2->refine_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 1];
1159: }
1160: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_z_hier_n) {
1161: dd2->coarsen_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 2];
1162: }
1163: dd2->refine_z_hier_n = dd->refine_z_hier_n;
1164: PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1165: PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);
1166: }
1167: if (dd->refine_y_hier) {
1168: if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1< dd->refine_y_hier_n) {
1169: dd2->refine_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 1];
1170: }
1171: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_y_hier_n) {
1172: dd2->coarsen_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 2];
1173: }
1174: dd2->refine_y_hier_n = dd->refine_y_hier_n;
1175: PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1176: PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);
1177: }
1178: if (dd->refine_x_hier) {
1179: if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_x_hier_n) {
1180: dd2->refine_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 1];
1181: }
1182: if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_x_hier_n) {
1183: dd2->coarsen_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 2];
1184: }
1185: dd2->refine_x_hier_n = dd->refine_x_hier_n;
1186: PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1187: PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);
1188: }
1190: /* copy vector type information */
1191: DMSetVecType(dmc2,dmf->vectype);
1193: dd2->lf = dd->lf;
1194: dd2->lj = dd->lj;
1196: dmc2->leveldown = dmf->leveldown + 1;
1197: dmc2->levelup = dmf->levelup;
1199: DMSetUp(dmc2);
1201: /* inject coordinates if they are set on the fine grid */
1202: if (dmf->coordinates) {
1203: DM cdaf,cdac;
1204: Vec coordsc,coordsf;
1205: Mat inject;
1206: VecScatter vscat;
1208: DMGetCoordinateDM(dmf,&cdaf);
1209: DMGetCoordinates(dmf,&coordsf);
1210: DMGetCoordinateDM(dmc2,&cdac);
1211: /* force creation of the coordinate vector */
1212: DMDASetUniformCoordinates(dmc2,0.0,1.0,0.0,1.0,0.0,1.0);
1213: DMGetCoordinates(dmc2,&coordsc);
1215: DMCreateInjection(cdac,cdaf,&inject);
1216: MatScatterGetVecScatter(inject,&vscat);
1217: VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1218: VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1219: MatDestroy(&inject);
1220: }
1222: for (i=0; i<dmf->bs; i++) {
1223: const char *fieldname;
1224: DMDAGetFieldName(dmf,i,&fieldname);
1225: DMDASetFieldName(dmc2,i,fieldname);
1226: }
1228: *dmc = dmc2;
1229: return(0);
1230: }
1232: PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1233: {
1235: PetscInt i,n,*refx,*refy,*refz;
1239: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1240: if (nlevels == 0) return(0);
1243: /* Get refinement factors, defaults taken from the coarse DMDA */
1244: PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1245: for (i=0; i<nlevels; i++) {
1246: DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1247: }
1248: n = nlevels;
1249: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1250: n = nlevels;
1251: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1252: n = nlevels;
1253: PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);
1255: DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1256: DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1257: for (i=1; i<nlevels; i++) {
1258: DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1259: DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1260: }
1261: PetscFree3(refx,refy,refz);
1262: return(0);
1263: }
1265: PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1266: {
1268: PetscInt i;
1272: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1273: if (nlevels == 0) return(0);
1275: DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1276: for (i=1; i<nlevels; i++) {
1277: DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1278: }
1279: return(0);
1280: }
1282: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm,PetscInt n,PetscReal *nodes)
1283: {
1285: PetscInt i,j,xs,xn,q;
1286: PetscScalar *xx;
1287: PetscReal h;
1288: Vec x;
1289: DM_DA *da = (DM_DA*)dm->data;
1292: if (da->bx != DM_BOUNDARY_PERIODIC) {
1293: DMDAGetInfo(dm,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1294: q = (q-1)/(n-1); /* number of spectral elements */
1295: h = 2.0/q;
1296: DMDAGetCorners(dm,&xs,NULL,NULL,&xn,NULL,NULL);
1297: xs = xs/(n-1);
1298: xn = xn/(n-1);
1299: DMDASetUniformCoordinates(dm,-1.,1.,0.,0.,0.,0.);
1300: DMGetCoordinates(dm,&x);
1301: DMDAVecGetArray(dm,x,&xx);
1303: /* loop over local spectral elements */
1304: for (j=xs; j<xs+xn; j++) {
1305: /*
1306: Except for the first process, each process starts on the second GLL point of the first element on that process
1307: */
1308: for (i= (j == xs && xs > 0)? 1 : 0; i<n; i++) {
1309: xx[j*(n-1) + i] = -1.0 + h*j + h*(nodes[i]+1.0)/2.;
1310: }
1311: }
1312: DMDAVecRestoreArray(dm,x,&xx);
1313: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for periodic");
1314: return(0);
1315: }
1317: /*@
1319: 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
1321: Collective on da
1323: Input Parameters:
1324: + da - the DMDA object
1325: - n - the number of GLL nodes
1326: - nodes - the GLL nodes
1328: Notes:
1329: the parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1330: on each process much be divisible by the number of GLL elements needed per process. This depends on whether the DM is
1331: periodic or not.
1333: Level: advanced
1335: .seealso: DMDACreate(), PetscDTGaussLobattoLegendreQuadrature(), DMGetCoordinates()
1336: @*/
1337: PetscErrorCode DMDASetGLLCoordinates(DM da,PetscInt n,PetscReal *nodes)
1338: {
1342: if (da->dim == 1) {
1343: DMDASetGLLCoordinates_1d(da,n,nodes);
1344: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for 2 or 3d");
1345: return(0);
1346: }
1348: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM da1,DM dm2,PetscBool *compatible,PetscBool *set)
1349: {
1351: DM_DA *dd1 = (DM_DA*)da1->data,*dd2;
1352: DM da2;
1353: DMType dmtype2;
1354: PetscBool isda,compatibleLocal;
1355: PetscInt i;
1358: if (!da1->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da1),PETSC_ERR_ARG_WRONGSTATE,"DMSetUp() must be called on first DM before DMGetCompatibility()");
1359: DMGetType(dm2,&dmtype2);
1360: PetscStrcmp(dmtype2,DMDA,&isda);
1361: if (isda) {
1362: da2 = dm2;
1363: dd2 = (DM_DA*)da2->data;
1364: if (!da2->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da2),PETSC_ERR_ARG_WRONGSTATE,"DMSetUp() must be called on second DM before DMGetCompatibility()");
1365: compatibleLocal = (PetscBool)(da1->dim == da2->dim);
1366: if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
1367: /* Global size ranks Boundary type */
1368: if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
1369: if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
1370: if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
1371: if (compatibleLocal) {
1372: for (i=0; i<dd1->m; ++i) {
1373: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i])); /* Local size */
1374: }
1375: }
1376: if (compatibleLocal && da1->dim > 1) {
1377: for (i=0; i<dd1->n; ++i) {
1378: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
1379: }
1380: }
1381: if (compatibleLocal && da1->dim > 2) {
1382: for (i=0; i<dd1->p; ++i) {
1383: compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
1384: }
1385: }
1386: *compatible = compatibleLocal;
1387: *set = PETSC_TRUE;
1388: } else {
1389: /* Decline to determine compatibility with other DM types */
1390: *set = PETSC_FALSE;
1391: }
1392: return(0);
1393: }