Actual source code: da.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/dmdaimpl.h>    /*I   "petscdmda.h"   I*/

  5: /*@
  6:   DMDASetSizes - Sets the global sizes

  8:   Logically Collective on DMDA

 10:   Input Parameters:
 11: + da - the DMDA
 12: . M - the global X size (or PETSC_DECIDE)
 13: . N - the global Y size (or PETSC_DECIDE)
 14: - P - the global Z size (or PETSC_DECIDE)

 16:   Level: intermediate

 18: .seealso: DMDAGetSize(), PetscSplitOwnership()
 19: @*/
 20: PetscErrorCode  DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
 21: {
 22:   DM_DA *dd = (DM_DA*)da->data;

 29:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");

 31:   dd->M = M;
 32:   dd->N = N;
 33:   dd->P = P;
 34:   return(0);
 35: }

 39: /*@
 40:   DMDASetNumProcs - Sets the number of processes in each dimension

 42:   Logically Collective on DMDA

 44:   Input Parameters:
 45: + da - the DMDA
 46: . m - the number of X procs (or PETSC_DECIDE)
 47: . n - the number of Y procs (or PETSC_DECIDE)
 48: - p - the number of Z procs (or PETSC_DECIDE)

 50:   Level: intermediate

 52: .seealso: DMDASetSizes(), DMDAGetSize(), PetscSplitOwnership()
 53: @*/
 54: PetscErrorCode  DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
 55: {
 56:   DM_DA          *dd = (DM_DA*)da->data;

 64:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
 65:   dd->m = m;
 66:   dd->n = n;
 67:   dd->p = p;
 68:   if (da->dim == 2) {
 69:     PetscMPIInt size;
 70:     MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
 71:     if ((dd->m > 0) && (dd->n < 0)) {
 72:       dd->n = size/dd->m;
 73:       if (dd->n*dd->m != size) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"%D processes in X direction not divisible into comm size %d",m,size);
 74:     }
 75:     if ((dd->n > 0) && (dd->m < 0)) {
 76:       dd->m = size/dd->n;
 77:       if (dd->n*dd->m != size) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"%D processes in Y direction not divisible into comm size %d",n,size);
 78:     }
 79:   }
 80:   return(0);
 81: }

 85: /*@
 86:   DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.

 88:   Not collective

 90:   Input Parameter:
 91: + da    - The DMDA
 92: - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC

 94:   Level: intermediate

 96: .keywords:  distributed array, periodicity
 97: .seealso: DMDACreate(), DMDestroy(), DMDA, DMBoundaryType
 98: @*/
 99: PetscErrorCode  DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)
100: {
101:   DM_DA *dd = (DM_DA*)da->data;

108:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
109:   dd->bx = bx;
110:   dd->by = by;
111:   dd->bz = bz;
112:   return(0);
113: }

117: /*@
118:   DMDASetDof - Sets the number of degrees of freedom per vertex

120:   Not collective

122:   Input Parameters:
123: + da  - The DMDA
124: - dof - Number of degrees of freedom

126:   Level: intermediate

128: .keywords:  distributed array, degrees of freedom
129: .seealso: DMDAGetDof(), DMDACreate(), DMDestroy(), DMDA
130: @*/
131: PetscErrorCode  DMDASetDof(DM da, PetscInt dof)
132: {
133:   DM_DA *dd = (DM_DA*)da->data;

138:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
139:   dd->w  = dof;
140:   da->bs = dof;
141:   return(0);
142: }

146: /*@
147:   DMDAGetDof - Gets the number of degrees of freedom per vertex

149:   Not collective

151:   Input Parameter:
152: . da  - The DMDA

154:   Output Parameter:
155: . dof - Number of degrees of freedom

157:   Level: intermediate

159: .keywords:  distributed array, degrees of freedom
160: .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
161: @*/
162: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
163: {
164:   DM_DA *dd = (DM_DA *) da->data;

169:   *dof = dd->w;
170:   return(0);
171: }

175: /*@
176:   DMDAGetOverlap - Gets the size of the per-processor overlap.

178:   Not collective

180:   Input Parameters:
181: . da  - The DMDA

183:   Output Parameters:
184: + x   - Overlap in the x direction
185: . y   - Overlap in the y direction
186: - z   - Overlap in the z direction

188:   Level: intermediate

190: .keywords:  distributed array, overlap, domain decomposition
191: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
192: @*/
193: PetscErrorCode  DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
194: {
195:   DM_DA *dd = (DM_DA*)da->data;

199:   if (x) *x = dd->xol;
200:   if (y) *y = dd->yol;
201:   if (z) *z = dd->zol;
202:   return(0);
203: }

207: /*@
208:   DMDASetOverlap - Sets the size of the per-processor overlap.

210:   Not collective

212:   Input Parameters:
213: + da  - The DMDA
214: . x   - Overlap in the x direction
215: . y   - Overlap in the y direction
216: - z   - Overlap in the z direction

218:   Level: intermediate

220: .keywords:  distributed array, overlap, domain decomposition
221: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
222: @*/
223: PetscErrorCode  DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
224: {
225:   DM_DA *dd = (DM_DA*)da->data;

232:   dd->xol = x;
233:   dd->yol = y;
234:   dd->zol = z;
235:   return(0);
236: }


241: /*@
242:   DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.

244:   Not collective

246:   Input Parameters:
247: . da  - The DMDA

249:   Output Parameters:
250: + Nsub   - Number of local subdomains created upon decomposition

252:   Level: intermediate

254: .keywords:  distributed array, domain decomposition
255: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
256: @*/
257: PetscErrorCode  DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
258: {
259:   DM_DA *dd = (DM_DA*)da->data;

263:   if (Nsub) *Nsub = dd->Nsub;
264:   return(0);
265: }

269: /*@
270:   DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.

272:   Not collective

274:   Input Parameters:
275: + da  - The DMDA
276: - Nsub - The number of local subdomains requested

278:   Level: intermediate

280: .keywords:  distributed array, domain decomposition
281: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
282: @*/
283: PetscErrorCode  DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
284: {
285:   DM_DA *dd = (DM_DA*)da->data;

290:   dd->Nsub = Nsub;
291:   return(0);
292: }

296: /*@
297:   DMDASetOffset - Sets the index offset of the DA.

299:   Collective on DA

301:   Input Parameter:
302: + da  - The DMDA
303: . xo  - The offset in the x direction
304: . yo  - The offset in the y direction
305: - zo  - The offset in the z direction

307:   Level: intermediate

309:   Notes: This is used primarily to overlap a computation on a local DA with that on a global DA without
310:   changing boundary conditions or subdomain features that depend upon the global offsets.

312: .keywords:  distributed array, degrees of freedom
313: .seealso: DMDAGetOffset(), DMDAVecGetArray()
314: @*/
315: PetscErrorCode  DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
316: {
318:   DM_DA          *dd = (DM_DA*)da->data;

328:   dd->xo = xo;
329:   dd->yo = yo;
330:   dd->zo = zo;
331:   dd->Mo = Mo;
332:   dd->No = No;
333:   dd->Po = Po;

335:   if (da->coordinateDM) {
336:     DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
337:   }
338:   return(0);
339: }

343: /*@
344:   DMDAGetOffset - Gets the index offset of the DA.

346:   Not collective

348:   Input Parameter:
349: . da  - The DMDA

351:   Output Parameters:
352: + xo  - The offset in the x direction
353: . yo  - The offset in the y direction
354: . zo  - The offset in the z direction
355: . Mo  - The global size in the x direction
356: . No  - The global size in the y direction
357: - Po  - The global size in the z direction

359:   Level: intermediate

361: .keywords:  distributed array, degrees of freedom
362: .seealso: DMDASetOffset(), DMDAVecGetArray()
363: @*/
364: PetscErrorCode  DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
365: {
366:   DM_DA *dd = (DM_DA*)da->data;

370:   if (xo) *xo = dd->xo;
371:   if (yo) *yo = dd->yo;
372:   if (zo) *zo = dd->zo;
373:   if (Mo) *Mo = dd->Mo;
374:   if (No) *No = dd->No;
375:   if (Po) *Po = dd->Po;
376:   return(0);
377: }

381: /*@
382:   DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.

384:   Not collective

386:   Input Parameter:
387: . da  - The DMDA

389:   Output Parameters:
390: + xs  - The start of the region in x
391: . ys  - The start of the region in y
392: . zs  - The start of the region in z
393: . xs  - The size of the region in x
394: . ys  - The size of the region in y
395: . zs  - The size of the region in z

397:   Level: intermediate

399: .keywords:  distributed array, degrees of freedom
400: .seealso: DMDAGetOffset(), DMDAVecGetArray()
401: @*/
402: PetscErrorCode  DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
403: {
404:   DM_DA          *dd = (DM_DA*)da->data;

408:   if (xs) *xs = dd->nonxs;
409:   if (ys) *ys = dd->nonys;
410:   if (zs) *zs = dd->nonzs;
411:   if (xm) *xm = dd->nonxm;
412:   if (ym) *ym = dd->nonym;
413:   if (zm) *zm = dd->nonzm;
414:   return(0);
415: }


420: /*@
421:   DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.

423:   Collective on DA

425:   Input Parameter:
426: + da  - The DMDA
427: . xs  - The start of the region in x
428: . ys  - The start of the region in y
429: . zs  - The start of the region in z
430: . xs  - The size of the region in x
431: . ys  - The size of the region in y
432: . zs  - The size of the region in z

434:   Level: intermediate

436: .keywords:  distributed array, degrees of freedom
437: .seealso: DMDAGetOffset(), DMDAVecGetArray()
438: @*/
439: PetscErrorCode  DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
440: {
441:   DM_DA          *dd = (DM_DA*)da->data;

451:   dd->nonxs = xs;
452:   dd->nonys = ys;
453:   dd->nonzs = zs;
454:   dd->nonxm = xm;
455:   dd->nonym = ym;
456:   dd->nonzm = zm;

458:   return(0);
459: }

463: /*@
464:   DMDASetStencilType - Sets the type of the communication stencil

466:   Logically Collective on DMDA

468:   Input Parameter:
469: + da    - The DMDA
470: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.

472:   Level: intermediate

474: .keywords:  distributed array, stencil
475: .seealso: DMDACreate(), DMDestroy(), DMDA
476: @*/
477: PetscErrorCode  DMDASetStencilType(DM da, DMDAStencilType stype)
478: {
479:   DM_DA *dd = (DM_DA*)da->data;

484:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
485:   dd->stencil_type = stype;
486:   return(0);
487: }

491: /*@
492:   DMDAGetStencilType - Gets the type of the communication stencil

494:   Not collective

496:   Input Parameter:
497: . da    - The DMDA

499:   Output Parameter:
500: . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.

502:   Level: intermediate

504: .keywords:  distributed array, stencil
505: .seealso: DMDACreate(), DMDestroy(), DMDA
506: @*/
507: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
508: {
509:   DM_DA *dd = (DM_DA*)da->data;

514:   *stype = dd->stencil_type;
515:   return(0);
516: }

520: /*@
521:   DMDASetStencilWidth - Sets the width of the communication stencil

523:   Logically Collective on DMDA

525:   Input Parameter:
526: + da    - The DMDA
527: - width - The stencil width

529:   Level: intermediate

531: .keywords:  distributed array, stencil
532: .seealso: DMDACreate(), DMDestroy(), DMDA
533: @*/
534: PetscErrorCode  DMDASetStencilWidth(DM da, PetscInt width)
535: {
536:   DM_DA *dd = (DM_DA*)da->data;

541:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
542:   dd->s = width;
543:   return(0);
544: }

548: /*@
549:   DMDAGetStencilWidth - Gets the width of the communication stencil

551:   Not collective

553:   Input Parameter:
554: . da    - The DMDA

556:   Output Parameter:
557: . width - The stencil width

559:   Level: intermediate

561: .keywords:  distributed array, stencil
562: .seealso: DMDACreate(), DMDestroy(), DMDA
563: @*/
564: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
565: {
566:   DM_DA *dd = (DM_DA *) da->data;

571:   *width = dd->s;
572:   return(0);
573: }

577: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
578: {
579:   PetscInt i,sum;

582:   if (M < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
583:   for (i=sum=0; i<m; i++) sum += lx[i];
584:   if (sum != M) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
585:   return(0);
586: }

590: /*@
591:   DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process

593:   Logically Collective on DMDA

595:   Input Parameter:
596: + da - The DMDA
597: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
598: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
599: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.

601:   Level: intermediate

603:   Note: these numbers are NOT multiplied by the number of dof per node.

605: .keywords:  distributed array
606: .seealso: DMDACreate(), DMDestroy(), DMDA
607: @*/
608: PetscErrorCode  DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
609: {
611:   DM_DA          *dd = (DM_DA*)da->data;

615:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
616:   if (lx) {
617:     if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
618:     DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
619:     if (!dd->lx) {
620:       PetscMalloc1(dd->m, &dd->lx);
621:     }
622:     PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
623:   }
624:   if (ly) {
625:     if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
626:     DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
627:     if (!dd->ly) {
628:       PetscMalloc1(dd->n, &dd->ly);
629:     }
630:     PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
631:   }
632:   if (lz) {
633:     if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
634:     DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
635:     if (!dd->lz) {
636:       PetscMalloc1(dd->p, &dd->lz);
637:     }
638:     PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
639:   }
640:   return(0);
641: }

645: /*@
646:        DMDASetInterpolationType - Sets the type of interpolation that will be
647:           returned by DMCreateInterpolation()

649:    Logically Collective on DMDA

651:    Input Parameter:
652: +  da - initial distributed array
653: .  ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms

655:    Level: intermediate

657:    Notes: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()

659: .keywords:  distributed array, interpolation

661: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
662: @*/
663: PetscErrorCode  DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
664: {
665:   DM_DA *dd = (DM_DA*)da->data;

670:   dd->interptype = ctype;
671:   return(0);
672: }

676: /*@
677:        DMDAGetInterpolationType - Gets the type of interpolation that will be
678:           used by DMCreateInterpolation()

680:    Not Collective

682:    Input Parameter:
683: .  da - distributed array

685:    Output Parameter:
686: .  ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)

688:    Level: intermediate

690: .keywords:  distributed array, interpolation

692: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
693: @*/
694: PetscErrorCode  DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
695: {
696:   DM_DA *dd = (DM_DA*)da->data;

701:   *ctype = dd->interptype;
702:   return(0);
703: }

707: /*@C
708:       DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
709:         processes neighbors.

711:     Not Collective

713:    Input Parameter:
714: .     da - the DMDA object

716:    Output Parameters:
717: .     ranks - the neighbors ranks, stored with the x index increasing most rapidly.
718:               this process itself is in the list

720:    Notes: In 2d the array is of length 9, in 3d of length 27
721:           Not supported in 1d
722:           Do not free the array, it is freed when the DMDA is destroyed.

724:    Fortran Notes: In fortran you must pass in an array of the appropriate length.

726:    Level: intermediate

728: @*/
729: PetscErrorCode  DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
730: {
731:   DM_DA *dd = (DM_DA*)da->data;

735:   *ranks = dd->neighbors;
736:   return(0);
737: }

741: /*@C
742:       DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process

744:     Not Collective

746:    Input Parameter:
747: .     da - the DMDA object

749:    Output Parameter:
750: +     lx - ownership along x direction (optional)
751: .     ly - ownership along y direction (optional)
752: -     lz - ownership along z direction (optional)

754:    Level: intermediate

756:     Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()

758:     In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
759:     eighth arguments from DMDAGetInfo()

761:      In C you should not free these arrays, nor change the values in them. They will only have valid values while the
762:     DMDA they came from still exists (has not been destroyed).

764:     These numbers are NOT multiplied by the number of dof per node.

766: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
767: @*/
768: PetscErrorCode  DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
769: {
770:   DM_DA *dd = (DM_DA*)da->data;

774:   if (lx) *lx = dd->lx;
775:   if (ly) *ly = dd->ly;
776:   if (lz) *lz = dd->lz;
777:   return(0);
778: }

782: /*@
783:      DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined

785:     Logically Collective on DMDA

787:   Input Parameters:
788: +    da - the DMDA object
789: .    refine_x - ratio of fine grid to coarse in x direction (2 by default)
790: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
791: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

793:   Options Database:
794: +  -da_refine_x - refinement ratio in x direction
795: .  -da_refine_y - refinement ratio in y direction
796: -  -da_refine_z - refinement ratio in z direction

798:   Level: intermediate

800:     Notes: Pass PETSC_IGNORE to leave a value unchanged

802: .seealso: DMRefine(), DMDAGetRefinementFactor()
803: @*/
804: PetscErrorCode  DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
805: {
806:   DM_DA *dd = (DM_DA*)da->data;


814:   if (refine_x > 0) dd->refine_x = refine_x;
815:   if (refine_y > 0) dd->refine_y = refine_y;
816:   if (refine_z > 0) dd->refine_z = refine_z;
817:   return(0);
818: }

822: /*@C
823:      DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined

825:     Not Collective

827:   Input Parameter:
828: .    da - the DMDA object

830:   Output Parameters:
831: +    refine_x - ratio of fine grid to coarse in x direction (2 by default)
832: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
833: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

835:   Level: intermediate

837:     Notes: Pass NULL for values you do not need

839: .seealso: DMRefine(), DMDASetRefinementFactor()
840: @*/
841: PetscErrorCode  DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
842: {
843:   DM_DA *dd = (DM_DA*)da->data;

847:   if (refine_x) *refine_x = dd->refine_x;
848:   if (refine_y) *refine_y = dd->refine_y;
849:   if (refine_z) *refine_z = dd->refine_z;
850:   return(0);
851: }

855: /*@C
856:      DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.

858:     Logically Collective on DMDA

860:   Input Parameters:
861: +    da - the DMDA object
862: -    f - the function that allocates the matrix for that specific DMDA

864:   Level: developer

866:    Notes: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
867:        the diagonal and off-diagonal blocks of the matrix

869: .seealso: DMCreateMatrix(), DMDASetBlockFills()
870: @*/
871: PetscErrorCode  DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
872: {
875:   da->ops->creatematrix = f;
876:   return(0);
877: }

881: /*
882:   Creates "balanced" ownership ranges after refinement, constrained by the need for the
883:   fine grid boundaries to fall within one stencil width of the coarse partition.

885:   Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
886: */
887: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
888: {
889:   PetscInt       i,totalc = 0,remaining,startc = 0,startf = 0;

893:   if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
894:   if (ratio == 1) {
895:     PetscMemcpy(lf,lc,m*sizeof(lc[0]));
896:     return(0);
897:   }
898:   for (i=0; i<m; i++) totalc += lc[i];
899:   remaining = (!periodic) + ratio * (totalc - (!periodic));
900:   for (i=0; i<m; i++) {
901:     PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
902:     if (i == m-1) lf[i] = want;
903:     else {
904:       const PetscInt nextc = startc + lc[i];
905:       /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
906:        * coarse stencil width of the first coarse node in the next subdomain. */
907:       while ((startf+want)/ratio < nextc - stencil_width) want++;
908:       /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
909:        * coarse stencil width of the last coarse node in the current subdomain. */
910:       while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
911:       /* Make sure all constraints are satisfied */
912:       if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
913:           || ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible refined ownership range");
914:     }
915:     lf[i]      = want;
916:     startc    += lc[i];
917:     startf    += lf[i];
918:     remaining -= lf[i];
919:   }
920:   return(0);
921: }

925: /*
926:   Creates "balanced" ownership ranges after coarsening, constrained by the need for the
927:   fine grid boundaries to fall within one stencil width of the coarse partition.

929:   Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
930: */
931: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
932: {
933:   PetscInt       i,totalf,remaining,startc,startf;

937:   if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
938:   if (ratio == 1) {
939:     PetscMemcpy(lc,lf,m*sizeof(lf[0]));
940:     return(0);
941:   }
942:   for (i=0,totalf=0; i<m; i++) totalf += lf[i];
943:   remaining = (!periodic) + (totalf - (!periodic)) / ratio;
944:   for (i=0,startc=0,startf=0; i<m; i++) {
945:     PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
946:     if (i == m-1) lc[i] = want;
947:     else {
948:       const PetscInt nextf = startf+lf[i];
949:       /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
950:        * node is within one stencil width. */
951:       while (nextf/ratio < startc+want-stencil_width) want--;
952:       /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
953:        * fine node is within one stencil width. */
954:       while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
955:       if (want < 0 || want > remaining
956:           || (nextf/ratio < startc+want-stencil_width) || ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible coarsened ownership range");
957:     }
958:     lc[i]      = want;
959:     startc    += lc[i];
960:     startf    += lf[i];
961:     remaining -= lc[i];
962:   }
963:   return(0);
964: }

968: PetscErrorCode  DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
969: {
971:   PetscInt       M,N,P,i,dim;
972:   DM             da2;
973:   DM_DA          *dd = (DM_DA*)da->data,*dd2;


979:   DMGetDimension(da, &dim);
980:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
981:     M = dd->refine_x*dd->M;
982:   } else {
983:     M = 1 + dd->refine_x*(dd->M - 1);
984:   }
985:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
986:     if (dim > 1) {
987:       N = dd->refine_y*dd->N;
988:     } else {
989:       N = 1;
990:     }
991:   } else {
992:     N = 1 + dd->refine_y*(dd->N - 1);
993:   }
994:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
995:     if (dim > 2) {
996:       P = dd->refine_z*dd->P;
997:     } else {
998:       P = 1;
999:     }
1000:   } else {
1001:     P = 1 + dd->refine_z*(dd->P - 1);
1002:   }
1003:   DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1004:   DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1005:   DMSetDimension(da2,dim);
1006:   DMDASetSizes(da2,M,N,P);
1007:   DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1008:   DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1009:   DMDASetDof(da2,dd->w);
1010:   DMDASetStencilType(da2,dd->stencil_type);
1011:   DMDASetStencilWidth(da2,dd->s);
1012:   if (dim == 3) {
1013:     PetscInt *lx,*ly,*lz;
1014:     PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1015:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
1016:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
1017:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
1018:     DMDASetOwnershipRanges(da2,lx,ly,lz);
1019:     PetscFree3(lx,ly,lz);
1020:   } else if (dim == 2) {
1021:     PetscInt *lx,*ly;
1022:     PetscMalloc2(dd->m,&lx,dd->n,&ly);
1023:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
1024:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
1025:     DMDASetOwnershipRanges(da2,lx,ly,NULL);
1026:     PetscFree2(lx,ly);
1027:   } else if (dim == 1) {
1028:     PetscInt *lx;
1029:     PetscMalloc1(dd->m,&lx);
1030:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
1031:     DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1032:     PetscFree(lx);
1033:   }
1034:   dd2 = (DM_DA*)da2->data;

1036:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
1037:   da2->ops->creatematrix = da->ops->creatematrix;
1038:   /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
1039:   da2->ops->getcoloring = da->ops->getcoloring;
1040:   dd2->interptype       = dd->interptype;

1042:   /* copy fill information if given */
1043:   if (dd->dfill) {
1044:     PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1045:     PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1046:   }
1047:   if (dd->ofill) {
1048:     PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1049:     PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1050:   }
1051:   /* copy the refine information */
1052:   dd2->coarsen_x = dd2->refine_x = dd->refine_x;
1053:   dd2->coarsen_y = dd2->refine_y = dd->refine_y;
1054:   dd2->coarsen_z = dd2->refine_z = dd->refine_z;

1056:   /* copy vector type information */
1057:   DMSetVecType(da2,da->vectype);

1059:   dd2->lf = dd->lf;
1060:   dd2->lj = dd->lj;

1062:   da2->leveldown = da->leveldown;
1063:   da2->levelup   = da->levelup + 1;

1065:   DMSetFromOptions(da2);
1066:   DMSetUp(da2);

1068:   /* interpolate coordinates if they are set on the coarse grid */
1069:   if (da->coordinates) {
1070:     DM  cdaf,cdac;
1071:     Vec coordsc,coordsf;
1072:     Mat II;

1074:     DMGetCoordinateDM(da,&cdac);
1075:     DMGetCoordinates(da,&coordsc);
1076:     DMGetCoordinateDM(da2,&cdaf);
1077:     /* force creation of the coordinate vector */
1078:     DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1079:     DMGetCoordinates(da2,&coordsf);
1080:     DMCreateInterpolation(cdac,cdaf,&II,NULL);
1081:     MatInterpolate(II,coordsc,coordsf);
1082:     MatDestroy(&II);
1083:   }

1085:   for (i=0; i<da->bs; i++) {
1086:     const char *fieldname;
1087:     DMDAGetFieldName(da,i,&fieldname);
1088:     DMDASetFieldName(da2,i,fieldname);
1089:   }

1091:   *daref = da2;
1092:   return(0);
1093: }


1098: PetscErrorCode  DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1099: {
1101:   PetscInt       M,N,P,i,dim;
1102:   DM             da2;
1103:   DM_DA          *dd = (DM_DA*)da->data,*dd2;


1109:   DMGetDimension(da, &dim);
1110:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1111:     M = dd->M / dd->coarsen_x;
1112:   } else {
1113:     M = 1 + (dd->M - 1) / dd->coarsen_x;
1114:   }
1115:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1116:     if (dim > 1) {
1117:       N = dd->N / dd->coarsen_y;
1118:     } else {
1119:       N = 1;
1120:     }
1121:   } else {
1122:     N = 1 + (dd->N - 1) / dd->coarsen_y;
1123:   }
1124:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1125:     if (dim > 2) {
1126:       P = dd->P / dd->coarsen_z;
1127:     } else {
1128:       P = 1;
1129:     }
1130:   } else {
1131:     P = 1 + (dd->P - 1) / dd->coarsen_z;
1132:   }
1133:   DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1134:   DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1135:   DMSetDimension(da2,dim);
1136:   DMDASetSizes(da2,M,N,P);
1137:   DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1138:   DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1139:   DMDASetDof(da2,dd->w);
1140:   DMDASetStencilType(da2,dd->stencil_type);
1141:   DMDASetStencilWidth(da2,dd->s);
1142:   if (dim == 3) {
1143:     PetscInt *lx,*ly,*lz;
1144:     PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1145:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1146:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1147:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1148:     DMDASetOwnershipRanges(da2,lx,ly,lz);
1149:     PetscFree3(lx,ly,lz);
1150:   } else if (dim == 2) {
1151:     PetscInt *lx,*ly;
1152:     PetscMalloc2(dd->m,&lx,dd->n,&ly);
1153:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1154:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1155:     DMDASetOwnershipRanges(da2,lx,ly,NULL);
1156:     PetscFree2(lx,ly);
1157:   } else if (dim == 1) {
1158:     PetscInt *lx;
1159:     PetscMalloc1(dd->m,&lx);
1160:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1161:     DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1162:     PetscFree(lx);
1163:   }
1164:   dd2 = (DM_DA*)da2->data;

1166:   /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1167:   /* da2->ops->createinterpolation = da->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1168:   da2->ops->creatematrix = da->ops->creatematrix;
1169:   da2->ops->getcoloring  = da->ops->getcoloring;
1170:   dd2->interptype        = dd->interptype;

1172:   /* copy fill information if given */
1173:   if (dd->dfill) {
1174:     PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1175:     PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1176:   }
1177:   if (dd->ofill) {
1178:     PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1179:     PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1180:   }
1181:   /* copy the refine information */
1182:   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1183:   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1184:   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;

1186:   /* copy vector type information */
1187:   DMSetVecType(da2,da->vectype);

1189:   dd2->lf = dd->lf;
1190:   dd2->lj = dd->lj;

1192:   da2->leveldown = da->leveldown + 1;
1193:   da2->levelup   = da->levelup;

1195:   DMSetFromOptions(da2);
1196:   DMSetUp(da2);

1198:   /* inject coordinates if they are set on the fine grid */
1199:   if (da->coordinates) {
1200:     DM         cdaf,cdac;
1201:     Vec        coordsc,coordsf;
1202:     Mat        inject;
1203:     VecScatter vscat;

1205:     DMGetCoordinateDM(da,&cdaf);
1206:     DMGetCoordinates(da,&coordsf);
1207:     DMGetCoordinateDM(da2,&cdac);
1208:     /* force creation of the coordinate vector */
1209:     DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1210:     DMGetCoordinates(da2,&coordsc);

1212:     DMCreateInjection(cdac,cdaf,&inject);
1213:     MatScatterGetVecScatter(inject,&vscat);
1214:     VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1215:     VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1216:     MatDestroy(&inject);
1217:   }

1219:   for (i=0; i<da->bs; i++) {
1220:     const char *fieldname;
1221:     DMDAGetFieldName(da,i,&fieldname);
1222:     DMDASetFieldName(da2,i,fieldname);
1223:   }

1225:   *daref = da2;
1226:   return(0);
1227: }

1231: PetscErrorCode  DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1232: {
1234:   PetscInt       i,n,*refx,*refy,*refz;

1238:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1239:   if (nlevels == 0) return(0);

1242:   /* Get refinement factors, defaults taken from the coarse DMDA */
1243:   PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1244:   for (i=0; i<nlevels; i++) {
1245:     DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1246:   }
1247:   n    = nlevels;
1248:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1249:   n    = nlevels;
1250:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1251:   n    = nlevels;
1252:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);

1254:   DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1255:   DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1256:   for (i=1; i<nlevels; i++) {
1257:     DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1258:     DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1259:   }
1260:   PetscFree3(refx,refy,refz);
1261:   return(0);
1262: }

1266: PetscErrorCode  DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1267: {
1269:   PetscInt       i;

1273:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1274:   if (nlevels == 0) return(0);
1276:   DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1277:   for (i=1; i<nlevels; i++) {
1278:     DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1279:   }
1280:   return(0);
1281: }