Actual source code: da.c
petsc-3.3-p7 2013-05-11
1: #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/
6: /*@
7: DMDASetDim - Sets the dimension
9: Collective on DMDA
11: Input Parameters:
12: + da - the DMDA
13: - dim - the dimension (or PETSC_DECIDE)
15: Level: intermediate
17: .seealso: DaGetDim(), DMDASetSizes()
18: @*/
19: PetscErrorCode DMDASetDim(DM da, PetscInt dim)
20: {
21: DM_DA *dd = (DM_DA*)da->data;
25: if (dd->dim > 0 && dim != dd->dim) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change DMDA dim from %D after it was set to %D",dd->dim,dim);
26: dd->dim = dim;
27: return(0);
28: }
32: /*@
33: DMDASetSizes - Sets the global sizes
35: Logically Collective on DMDA
37: Input Parameters:
38: + da - the DMDA
39: . M - the global X size (or PETSC_DECIDE)
40: . N - the global Y size (or PETSC_DECIDE)
41: - P - the global Z size (or PETSC_DECIDE)
43: Level: intermediate
45: .seealso: DMDAGetSize(), PetscSplitOwnership()
46: @*/
47: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
48: {
49: DM_DA *dd = (DM_DA*)da->data;
56: if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
58: dd->M = M;
59: dd->N = N;
60: dd->P = P;
61: return(0);
62: }
66: /*@
67: DMDASetNumProcs - Sets the number of processes in each dimension
69: Logically Collective on DMDA
71: Input Parameters:
72: + da - the DMDA
73: . m - the number of X procs (or PETSC_DECIDE)
74: . n - the number of Y procs (or PETSC_DECIDE)
75: - p - the number of Z procs (or PETSC_DECIDE)
77: Level: intermediate
79: .seealso: DMDASetSizes(), DMDAGetSize(), PetscSplitOwnership()
80: @*/
81: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
82: {
83: DM_DA *dd = (DM_DA*)da->data;
90: if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
91: dd->m = m;
92: dd->n = n;
93: dd->p = p;
94: return(0);
95: }
99: /*@
100: DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.
102: Not collective
104: Input Parameter:
105: + da - The DMDA
106: - bx,by,bz - One of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC
108: Level: intermediate
110: .keywords: distributed array, periodicity
111: .seealso: DMDACreate(), DMDestroy(), DMDA, DMDABoundaryType
112: @*/
113: PetscErrorCode DMDASetBoundaryType(DM da,DMDABoundaryType bx,DMDABoundaryType by,DMDABoundaryType bz)
114: {
115: DM_DA *dd = (DM_DA*)da->data;
122: if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
123: dd->bx = bx;
124: dd->by = by;
125: dd->bz = bz;
126: return(0);
127: }
131: /*@
132: DMDASetDof - Sets the number of degrees of freedom per vertex
134: Not collective
136: Input Parameter:
137: + da - The DMDA
138: - dof - Number of degrees of freedom
140: Level: intermediate
142: .keywords: distributed array, degrees of freedom
143: .seealso: DMDACreate(), DMDestroy(), DMDA
144: @*/
145: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
146: {
147: DM_DA *dd = (DM_DA*)da->data;
152: if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
153: dd->w = dof;
154: da->bs = dof;
155: return(0);
156: }
160: /*@
161: DMDASetStencilType - Sets the type of the communication stencil
163: Logically Collective on DMDA
165: Input Parameter:
166: + da - The DMDA
167: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
169: Level: intermediate
171: .keywords: distributed array, stencil
172: .seealso: DMDACreate(), DMDestroy(), DMDA
173: @*/
174: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
175: {
176: DM_DA *dd = (DM_DA*)da->data;
181: if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
182: dd->stencil_type = stype;
183: return(0);
184: }
188: /*@
189: DMDASetStencilWidth - Sets the width of the communication stencil
191: Logically Collective on DMDA
193: Input Parameter:
194: + da - The DMDA
195: - width - The stencil width
197: Level: intermediate
199: .keywords: distributed array, stencil
200: .seealso: DMDACreate(), DMDestroy(), DMDA
201: @*/
202: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
203: {
204: DM_DA *dd = (DM_DA*)da->data;
209: if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
210: dd->s = width;
211: return(0);
212: }
216: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
217: {
218: PetscInt i,sum;
221: if (M < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
222: for (i=sum=0; i<m; i++) sum += lx[i];
223: if (sum != M) SETERRQ2(((PetscObject)da)->comm,PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
224: return(0);
225: }
229: /*@
230: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
232: Logically Collective on DMDA
234: Input Parameter:
235: + da - The DMDA
236: . lx - array containing number of nodes in the X direction on each process, or PETSC_NULL. If non-null, must be of length da->m
237: . ly - array containing number of nodes in the Y direction on each process, or PETSC_NULL. If non-null, must be of length da->n
238: - lz - array containing number of nodes in the Z direction on each process, or PETSC_NULL. If non-null, must be of length da->p.
240: Level: intermediate
242: .keywords: distributed array
243: .seealso: DMDACreate(), DMDestroy(), DMDA
244: @*/
245: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
246: {
248: DM_DA *dd = (DM_DA*)da->data;
249:
252: if (da->setupcalled) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
253: if (lx) {
254: if (dd->m < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
255: DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
256: if (!dd->lx) {
257: PetscMalloc(dd->m*sizeof(PetscInt), &dd->lx);
258: }
259: PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
260: }
261: if (ly) {
262: if (dd->n < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
263: DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
264: if (!dd->ly) {
265: PetscMalloc(dd->n*sizeof(PetscInt), &dd->ly);
266: }
267: PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
268: }
269: if (lz) {
270: if (dd->p < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
271: DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
272: if (!dd->lz) {
273: PetscMalloc(dd->p*sizeof(PetscInt), &dd->lz);
274: }
275: PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
276: }
277: return(0);
278: }
282: /*@
283: DMDASetInterpolationType - Sets the type of interpolation that will be
284: returned by DMCreateInterpolation()
286: Logically Collective on DMDA
288: Input Parameter:
289: + da - initial distributed array
290: . ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
292: Level: intermediate
294: Notes: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
296: .keywords: distributed array, interpolation
298: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
299: @*/
300: PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
301: {
302: DM_DA *dd = (DM_DA*)da->data;
307: dd->interptype = ctype;
308: return(0);
309: }
313: /*@
314: DMDAGetInterpolationType - Gets the type of interpolation that will be
315: used by DMCreateInterpolation()
317: Not Collective
319: Input Parameter:
320: . da - distributed array
322: Output Parameter:
323: . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
325: Level: intermediate
327: .keywords: distributed array, interpolation
329: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
330: @*/
331: PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
332: {
333: DM_DA *dd = (DM_DA*)da->data;
338: *ctype = dd->interptype;
339: return(0);
340: }
344: /*@C
345: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
346: processes neighbors.
348: Not Collective
350: Input Parameter:
351: . da - the DMDA object
353: Output Parameters:
354: . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
355: this process itself is in the list
357: Notes: In 2d the array is of length 9, in 3d of length 27
358: Not supported in 1d
359: Do not free the array, it is freed when the DMDA is destroyed.
361: Fortran Notes: In fortran you must pass in an array of the appropriate length.
363: Level: intermediate
365: @*/
366: PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
367: {
368: DM_DA *dd = (DM_DA*)da->data;
371: *ranks = dd->neighbors;
372: return(0);
373: }
377: /*@C
378: DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
380: Not Collective
382: Input Parameter:
383: . da - the DMDA object
385: Output Parameter:
386: + lx - ownership along x direction (optional)
387: . ly - ownership along y direction (optional)
388: - lz - ownership along z direction (optional)
390: Level: intermediate
392: Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
394: In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
395: eighth arguments from DMDAGetInfo()
397: In C you should not free these arrays, nor change the values in them. They will only have valid values while the
398: DMDA they came from still exists (has not been destroyed).
400: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
401: @*/
402: PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
403: {
404: DM_DA *dd = (DM_DA*)da->data;
408: if (lx) *lx = dd->lx;
409: if (ly) *ly = dd->ly;
410: if (lz) *lz = dd->lz;
411: return(0);
412: }
416: /*@
417: DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
419: Logically Collective on DMDA
421: Input Parameters:
422: + da - the DMDA object
423: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
424: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
425: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
427: Options Database:
428: + -da_refine_x - refinement ratio in x direction
429: . -da_refine_y - refinement ratio in y direction
430: - -da_refine_z - refinement ratio in z direction
432: Level: intermediate
434: Notes: Pass PETSC_IGNORE to leave a value unchanged
436: .seealso: DMRefine(), DMDAGetRefinementFactor()
437: @*/
438: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
439: {
440: DM_DA *dd = (DM_DA*)da->data;
448: if (refine_x > 0) dd->refine_x = refine_x;
449: if (refine_y > 0) dd->refine_y = refine_y;
450: if (refine_z > 0) dd->refine_z = refine_z;
451: return(0);
452: }
456: /*@C
457: DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
459: Not Collective
461: Input Parameter:
462: . da - the DMDA object
464: Output Parameters:
465: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
466: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
467: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
469: Level: intermediate
471: Notes: Pass PETSC_NULL for values you do not need
473: .seealso: DMRefine(), DMDASetRefinementFactor()
474: @*/
475: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
476: {
477: DM_DA *dd = (DM_DA*)da->data;
481: if (refine_x) *refine_x = dd->refine_x;
482: if (refine_y) *refine_y = dd->refine_y;
483: if (refine_z) *refine_z = dd->refine_z;
484: return(0);
485: }
489: /*@C
490: DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
492: Logically Collective on DMDA
494: Input Parameters:
495: + da - the DMDA object
496: - f - the function that allocates the matrix for that specific DMDA
498: Level: developer
500: Notes: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
501: the diagonal and off-diagonal blocks of the matrix
503: .seealso: DMCreateMatrix(), DMDASetBlockFills()
504: @*/
505: PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, const MatType,Mat*))
506: {
509: da->ops->creatematrix = f;
510: return(0);
511: }
515: /*
516: Creates "balanced" ownership ranges after refinement, constrained by the need for the
517: fine grid boundaries to fall within one stencil width of the coarse partition.
519: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
520: */
521: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
522: {
523: PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
527: if (ratio < 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
528: if (ratio == 1) {
529: PetscMemcpy(lf,lc,m*sizeof(lc[0]));
530: return(0);
531: }
532: for (i=0; i<m; i++) totalc += lc[i];
533: remaining = (!periodic) + ratio * (totalc - (!periodic));
534: for (i=0; i<m; i++) {
535: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
536: if (i == m-1) lf[i] = want;
537: else {
538: const PetscInt nextc = startc + lc[i];
539: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
540: * coarse stencil width of the first coarse node in the next subdomain. */
541: while ((startf+want)/ratio < nextc - stencil_width) want++;
542: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
543: * coarse stencil width of the last coarse node in the current subdomain. */
544: while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
545: /* Make sure all constraints are satisfied */
546: if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
547: || ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Could not find a compatible refined ownership range");
548: }
549: lf[i] = want;
550: startc += lc[i];
551: startf += lf[i];
552: remaining -= lf[i];
553: }
554: return(0);
555: }
559: /*
560: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
561: fine grid boundaries to fall within one stencil width of the coarse partition.
563: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
564: */
565: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
566: {
567: PetscInt i,totalf,remaining,startc,startf;
571: if (ratio < 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
572: if (ratio == 1) {
573: PetscMemcpy(lc,lf,m*sizeof(lf[0]));
574: return(0);
575: }
576: for (i=0,totalf=0; i<m; i++) totalf += lf[i];
577: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
578: for (i=0,startc=0,startf=0; i<m; i++) {
579: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
580: if (i == m-1) lc[i] = want;
581: else {
582: const PetscInt nextf = startf+lf[i];
583: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
584: * node is within one stencil width. */
585: while (nextf/ratio < startc+want-stencil_width) want--;
586: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
587: * fine node is within one stencil width. */
588: while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
589: if (want < 0 || want > remaining
590: || (nextf/ratio < startc+want-stencil_width) || ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width)) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Could not find a compatible coarsened ownership range");
591: }
592: lc[i] = want;
593: startc += lc[i];
594: startf += lf[i];
595: remaining -= lc[i];
596: }
597: return(0);
598: }
602: PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
603: {
605: PetscInt M,N,P,i;
606: DM da2;
607: DM_DA *dd = (DM_DA*)da->data,*dd2;
613: if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
614: M = dd->refine_x*dd->M;
615: } else {
616: M = 1 + dd->refine_x*(dd->M - 1);
617: }
618: if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
619: if (dd->dim > 1) {
620: N = dd->refine_y*dd->N;
621: } else {
622: N = 1;
623: }
624: } else {
625: N = 1 + dd->refine_y*(dd->N - 1);
626: }
627: if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
628: if (dd->dim > 2) {
629: P = dd->refine_z*dd->P;
630: } else {
631: P = 1;
632: }
633: } else {
634: P = 1 + dd->refine_z*(dd->P - 1);
635: }
636: DMDACreate(((PetscObject)da)->comm,&da2);
637: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
638: DMDASetDim(da2,dd->dim);
639: DMDASetSizes(da2,M,N,P);
640: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
641: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
642: DMDASetDof(da2,dd->w);
643: DMDASetStencilType(da2,dd->stencil_type);
644: DMDASetStencilWidth(da2,dd->s);
645: if (dd->dim == 3) {
646: PetscInt *lx,*ly,*lz;
647: PetscMalloc3(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly,dd->p,PetscInt,&lz);
648: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
649: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
650: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
651: DMDASetOwnershipRanges(da2,lx,ly,lz);
652: PetscFree3(lx,ly,lz);
653: } else if (dd->dim == 2) {
654: PetscInt *lx,*ly;
655: PetscMalloc2(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly);
656: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
657: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
658: DMDASetOwnershipRanges(da2,lx,ly,PETSC_NULL);
659: PetscFree2(lx,ly);
660: } else if (dd->dim == 1) {
661: PetscInt *lx;
662: PetscMalloc(dd->m*sizeof(PetscInt),&lx);
663: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
664: DMDASetOwnershipRanges(da2,lx,PETSC_NULL,PETSC_NULL);
665: PetscFree(lx);
666: }
667: dd2 = (DM_DA*)da2->data;
669: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
670: da2->ops->creatematrix = da->ops->creatematrix;
671: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
672: da2->ops->getcoloring = da->ops->getcoloring;
673: dd2->interptype = dd->interptype;
674:
675: /* copy fill information if given */
676: if (dd->dfill) {
677: PetscMalloc((dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->dfill);
678: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
679: }
680: if (dd->ofill) {
681: PetscMalloc((dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->ofill);
682: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
683: }
684: /* copy the refine information */
685: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
686: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
687: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
689: /* copy vector type information */
690: PetscFree(da2->vectype);
691: PetscStrallocpy(da->vectype,&da2->vectype);
693: dd2->lf = dd->lf;
694: dd2->lj = dd->lj;
696: da2->leveldown = da->leveldown;
697: da2->levelup = da->levelup + 1;
698: DMSetFromOptions(da2);
699: DMSetUp(da2);
700: DMView_DA_Private(da2);
702: /* interpolate coordinates if they are set on the coarse grid */
703: if (dd->coordinates) {
704: DM cdaf,cdac;
705: Vec coordsc,coordsf;
706: Mat II;
707:
708: DMDAGetCoordinateDA(da,&cdac);
709: DMDAGetCoordinates(da,&coordsc);
710: DMDAGetCoordinateDA(da2,&cdaf);
711: /* force creation of the coordinate vector */
712: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
713: DMDAGetCoordinates(da2,&coordsf);
714: DMCreateInterpolation(cdac,cdaf,&II,PETSC_NULL);
715: MatInterpolate(II,coordsc,coordsf);
716: MatDestroy(&II);
717: }
719: for (i=0; i<da->bs; i++) {
720: const char *fieldname;
721: DMDAGetFieldName(da,i,&fieldname);
722: DMDASetFieldName(da2,i,fieldname);
723: }
725: *daref = da2;
726: return(0);
727: }
732: PetscErrorCode DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
733: {
735: PetscInt M,N,P,i;
736: DM da2;
737: DM_DA *dd = (DM_DA*)da->data,*dd2;
743: if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
744: M = dd->M / dd->coarsen_x;
745: } else {
746: M = 1 + (dd->M - 1) / dd->coarsen_x;
747: }
748: if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
749: if (dd->dim > 1) {
750: N = dd->N / dd->coarsen_y;
751: } else {
752: N = 1;
753: }
754: } else {
755: N = 1 + (dd->N - 1) / dd->coarsen_y;
756: }
757: if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0){
758: if (dd->dim > 2) {
759: P = dd->P / dd->coarsen_z;
760: } else {
761: P = 1;
762: }
763: } else {
764: P = 1 + (dd->P - 1) / dd->coarsen_z;
765: }
766: DMDACreate(((PetscObject)da)->comm,&da2);
767: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
768: DMDASetDim(da2,dd->dim);
769: DMDASetSizes(da2,M,N,P);
770: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
771: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
772: DMDASetDof(da2,dd->w);
773: DMDASetStencilType(da2,dd->stencil_type);
774: DMDASetStencilWidth(da2,dd->s);
775: if (dd->dim == 3) {
776: PetscInt *lx,*ly,*lz;
777: PetscMalloc3(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly,dd->p,PetscInt,&lz);
778: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
779: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
780: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
781: DMDASetOwnershipRanges(da2,lx,ly,lz);
782: PetscFree3(lx,ly,lz);
783: } else if (dd->dim == 2) {
784: PetscInt *lx,*ly;
785: PetscMalloc2(dd->m,PetscInt,&lx,dd->n,PetscInt,&ly);
786: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
787: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
788: DMDASetOwnershipRanges(da2,lx,ly,PETSC_NULL);
789: PetscFree2(lx,ly);
790: } else if (dd->dim == 1) {
791: PetscInt *lx;
792: PetscMalloc(dd->m*sizeof(PetscInt),&lx);
793: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
794: DMDASetOwnershipRanges(da2,lx,PETSC_NULL,PETSC_NULL);
795: PetscFree(lx);
796: }
797: dd2 = (DM_DA*)da2->data;
799: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
800: /* da2->ops->createinterpolation = da->ops->createinterpolation; copying this one causes trouble for DMSetVI */
801: da2->ops->creatematrix = da->ops->creatematrix;
802: da2->ops->getcoloring = da->ops->getcoloring;
803: dd2->interptype = dd->interptype;
804:
805: /* copy fill information if given */
806: if (dd->dfill) {
807: PetscMalloc((dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->dfill);
808: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
809: }
810: if (dd->ofill) {
811: PetscMalloc((dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt),&dd2->ofill);
812: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
813: }
814: /* copy the refine information */
815: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
816: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
817: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
819: /* copy vector type information */
820: PetscFree(da2->vectype);
821: PetscStrallocpy(da->vectype,&da2->vectype);
823: dd2->lf = dd->lf;
824: dd2->lj = dd->lj;
826: da2->leveldown = da->leveldown + 1;
827: da2->levelup = da->levelup;
828: DMSetFromOptions(da2);
829: DMSetUp(da2);
830: DMView_DA_Private(da2);
832: /* inject coordinates if they are set on the fine grid */
833: if (dd->coordinates) {
834: DM cdaf,cdac;
835: Vec coordsc,coordsf;
836: VecScatter inject;
837:
838: DMDAGetCoordinateDA(da,&cdaf);
839: DMDAGetCoordinates(da,&coordsf);
840: DMDAGetCoordinateDA(da2,&cdac);
841: /* force creation of the coordinate vector */
842: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
843: DMDAGetCoordinates(da2,&coordsc);
844:
845: DMCreateInjection(cdac,cdaf,&inject);
846: VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
847: VecScatterEnd(inject ,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
848: VecScatterDestroy(&inject);
849: }
851: for (i=0; i<da->bs; i++) {
852: const char *fieldname;
853: DMDAGetFieldName(da,i,&fieldname);
854: DMDASetFieldName(da2,i,fieldname);
855: }
857: *daref = da2;
858: return(0);
859: }
863: PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
864: {
866: PetscInt i,n,*refx,*refy,*refz;
870: if (nlevels < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
871: if (nlevels == 0) return(0);
874: /* Get refinement factors, defaults taken from the coarse DMDA */
875: PetscMalloc3(nlevels,PetscInt,&refx,nlevels,PetscInt,&refy,nlevels,PetscInt,&refz);
876: for (i=0; i<nlevels; i++) {
877: DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
878: }
879: n = nlevels;
880: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,PETSC_NULL);
881: n = nlevels;
882: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,PETSC_NULL);
883: n = nlevels;
884: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,PETSC_NULL);
886: DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
887: DMRefine(da,((PetscObject)da)->comm,&daf[0]);
888: for (i=1; i<nlevels; i++) {
889: DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
890: DMRefine(daf[i-1],((PetscObject)da)->comm,&daf[i]);
891: }
892: PetscFree3(refx,refy,refz);
893: return(0);
894: }
898: PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
899: {
901: PetscInt i;
905: if (nlevels < 0) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
906: if (nlevels == 0) return(0);
908: DMCoarsen(da,((PetscObject)da)->comm,&dac[0]);
909: for (i=1; i<nlevels; i++) {
910: DMCoarsen(dac[i-1],((PetscObject)da)->comm,&dac[i]);
911: }
912: return(0);
913: }