Actual source code: da.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  1: #include <petsc-private/dmdaimpl.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(PetscObjectComm((PetscObject)da),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(PetscObjectComm((PetscObject)da),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;

 91:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
 92:   dd->m = m;
 93:   dd->n = n;
 94:   dd->p = p;
 95:   if (dd->dim == 2) {
 96:     PetscMPIInt size;
 97:     MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
 98:     if ((dd->m > 0) && (dd->n < 0)) {
 99:       dd->n = size/dd->m;
100:       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);
101:     }
102:     if ((dd->n > 0) && (dd->m < 0)) {
103:       dd->m = size/dd->n;
104:       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);
105:     }
106:   }
107:   return(0);
108: }

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

115:   Not collective

117:   Input Parameter:
118: + da    - The DMDA
119: - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC

121:   Level: intermediate

123: .keywords:  distributed array, periodicity
124: .seealso: DMDACreate(), DMDestroy(), DMDA, DMBoundaryType
125: @*/
126: PetscErrorCode  DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)
127: {
128:   DM_DA *dd = (DM_DA*)da->data;

135:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
136:   dd->bx = bx;
137:   dd->by = by;
138:   dd->bz = bz;
139:   return(0);
140: }

144: /*@
145:   DMDASetDof - Sets the number of degrees of freedom per vertex

147:   Not collective

149:   Input Parameter:
150: + da  - The DMDA
151: - dof - Number of degrees of freedom

153:   Level: intermediate

155: .keywords:  distributed array, degrees of freedom
156: .seealso: DMDACreate(), DMDestroy(), DMDA
157: @*/
158: PetscErrorCode  DMDASetDof(DM da, PetscInt dof)
159: {
160:   DM_DA *dd = (DM_DA*)da->data;

165:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
166:   dd->w  = dof;
167:   da->bs = dof;
168:   return(0);
169: }

173: /*@
174:   DMDAGetOverlap - Gets the size of the per-processor overlap.

176:   Not collective

178:   Input Parameters:
179: . da  - The DMDA

181:   Output Parameters:
182: + x   - Overlap in the x direction
183: . y   - Overlap in the y direction
184: - z   - Overlap in the z direction

186:   Level: intermediate

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

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

205: /*@
206:   DMDASetOverlap - Sets the size of the per-processor overlap.

208:   Not collective

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

216:   Level: intermediate

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

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


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

242:   Not collective

244:   Input Parameters:
245: . da  - The DMDA

247:   Output Parameters:
248: + Nsub   - Number of local subdomains created upon decomposition

250:   Level: intermediate

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

261:   if (Nsub) *Nsub = dd->Nsub;
262:   return(0);
263: }

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

270:   Not collective

272:   Input Parameters:
273: + da  - The DMDA
274: - Nsub - The number of local subdomains requested

276:   Level: intermediate

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

288:   dd->Nsub = Nsub;
289:   return(0);
290: }

294: /*@
295:   DMDASetOffset - Sets the index offset of the DA.

297:   Collective on DA

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

305:   Level: intermediate

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

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

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

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

341: /*@
342:   DMDAGetOffset - Gets the index offset of the DA.

344:   Not collective

346:   Input Parameter:
347: . da  - The DMDA

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

357:   Level: intermediate

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

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

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

382:   Not collective

384:   Input Parameter:
385: . da  - The DMDA

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

395:   Level: intermediate

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

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


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

421:   Collective on DA

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

432:   Level: intermediate

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

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

456:   return(0);
457: }

461: /*@
462:   DMDASetStencilType - Sets the type of the communication stencil

464:   Logically Collective on DMDA

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

470:   Level: intermediate

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

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

489: /*@
490:   DMDASetStencilWidth - Sets the width of the communication stencil

492:   Logically Collective on DMDA

494:   Input Parameter:
495: + da    - The DMDA
496: - width - The stencil width

498:   Level: intermediate

500: .keywords:  distributed array, stencil
501: .seealso: DMDACreate(), DMDestroy(), DMDA
502: @*/
503: PetscErrorCode  DMDASetStencilWidth(DM da, PetscInt width)
504: {
505:   DM_DA *dd = (DM_DA*)da->data;

510:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
511:   dd->s = width;
512:   return(0);
513: }

517: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
518: {
519:   PetscInt i,sum;

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

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

533:   Logically Collective on DMDA

535:   Input Parameter:
536: + da - The DMDA
537: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
538: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
539: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.

541:   Level: intermediate

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

545: .keywords:  distributed array
546: .seealso: DMDACreate(), DMDestroy(), DMDA
547: @*/
548: PetscErrorCode  DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
549: {
551:   DM_DA          *dd = (DM_DA*)da->data;

555:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
556:   if (lx) {
557:     if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
558:     DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
559:     if (!dd->lx) {
560:       PetscMalloc1(dd->m, &dd->lx);
561:     }
562:     PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
563:   }
564:   if (ly) {
565:     if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
566:     DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
567:     if (!dd->ly) {
568:       PetscMalloc1(dd->n, &dd->ly);
569:     }
570:     PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
571:   }
572:   if (lz) {
573:     if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
574:     DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
575:     if (!dd->lz) {
576:       PetscMalloc1(dd->p, &dd->lz);
577:     }
578:     PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
579:   }
580:   return(0);
581: }

585: /*@
586:        DMDASetInterpolationType - Sets the type of interpolation that will be
587:           returned by DMCreateInterpolation()

589:    Logically Collective on DMDA

591:    Input Parameter:
592: +  da - initial distributed array
593: .  ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms

595:    Level: intermediate

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

599: .keywords:  distributed array, interpolation

601: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
602: @*/
603: PetscErrorCode  DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
604: {
605:   DM_DA *dd = (DM_DA*)da->data;

610:   dd->interptype = ctype;
611:   return(0);
612: }

616: /*@
617:        DMDAGetInterpolationType - Gets the type of interpolation that will be
618:           used by DMCreateInterpolation()

620:    Not Collective

622:    Input Parameter:
623: .  da - distributed array

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

628:    Level: intermediate

630: .keywords:  distributed array, interpolation

632: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
633: @*/
634: PetscErrorCode  DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
635: {
636:   DM_DA *dd = (DM_DA*)da->data;

641:   *ctype = dd->interptype;
642:   return(0);
643: }

647: /*@C
648:       DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
649:         processes neighbors.

651:     Not Collective

653:    Input Parameter:
654: .     da - the DMDA object

656:    Output Parameters:
657: .     ranks - the neighbors ranks, stored with the x index increasing most rapidly.
658:               this process itself is in the list

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

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

666:    Level: intermediate

668: @*/
669: PetscErrorCode  DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
670: {
671:   DM_DA *dd = (DM_DA*)da->data;

675:   *ranks = dd->neighbors;
676:   return(0);
677: }

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

684:     Not Collective

686:    Input Parameter:
687: .     da - the DMDA object

689:    Output Parameter:
690: +     lx - ownership along x direction (optional)
691: .     ly - ownership along y direction (optional)
692: -     lz - ownership along z direction (optional)

694:    Level: intermediate

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

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

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

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

706: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
707: @*/
708: PetscErrorCode  DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
709: {
710:   DM_DA *dd = (DM_DA*)da->data;

714:   if (lx) *lx = dd->lx;
715:   if (ly) *ly = dd->ly;
716:   if (lz) *lz = dd->lz;
717:   return(0);
718: }

722: /*@
723:      DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined

725:     Logically Collective on DMDA

727:   Input Parameters:
728: +    da - the DMDA object
729: .    refine_x - ratio of fine grid to coarse in x direction (2 by default)
730: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
731: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

733:   Options Database:
734: +  -da_refine_x - refinement ratio in x direction
735: .  -da_refine_y - refinement ratio in y direction
736: -  -da_refine_z - refinement ratio in z direction

738:   Level: intermediate

740:     Notes: Pass PETSC_IGNORE to leave a value unchanged

742: .seealso: DMRefine(), DMDAGetRefinementFactor()
743: @*/
744: PetscErrorCode  DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
745: {
746:   DM_DA *dd = (DM_DA*)da->data;


754:   if (refine_x > 0) dd->refine_x = refine_x;
755:   if (refine_y > 0) dd->refine_y = refine_y;
756:   if (refine_z > 0) dd->refine_z = refine_z;
757:   return(0);
758: }

762: /*@C
763:      DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined

765:     Not Collective

767:   Input Parameter:
768: .    da - the DMDA object

770:   Output Parameters:
771: +    refine_x - ratio of fine grid to coarse in x direction (2 by default)
772: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
773: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

775:   Level: intermediate

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

779: .seealso: DMRefine(), DMDASetRefinementFactor()
780: @*/
781: PetscErrorCode  DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
782: {
783:   DM_DA *dd = (DM_DA*)da->data;

787:   if (refine_x) *refine_x = dd->refine_x;
788:   if (refine_y) *refine_y = dd->refine_y;
789:   if (refine_z) *refine_z = dd->refine_z;
790:   return(0);
791: }

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

798:     Logically Collective on DMDA

800:   Input Parameters:
801: +    da - the DMDA object
802: -    f - the function that allocates the matrix for that specific DMDA

804:   Level: developer

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

809: .seealso: DMCreateMatrix(), DMDASetBlockFills()
810: @*/
811: PetscErrorCode  DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
812: {
815:   da->ops->creatematrix = f;
816:   return(0);
817: }

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

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

833:   if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
834:   if (ratio == 1) {
835:     PetscMemcpy(lf,lc,m*sizeof(lc[0]));
836:     return(0);
837:   }
838:   for (i=0; i<m; i++) totalc += lc[i];
839:   remaining = (!periodic) + ratio * (totalc - (!periodic));
840:   for (i=0; i<m; i++) {
841:     PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
842:     if (i == m-1) lf[i] = want;
843:     else {
844:       const PetscInt nextc = startc + lc[i];
845:       /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
846:        * coarse stencil width of the first coarse node in the next subdomain. */
847:       while ((startf+want)/ratio < nextc - stencil_width) want++;
848:       /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
849:        * coarse stencil width of the last coarse node in the current subdomain. */
850:       while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
851:       /* Make sure all constraints are satisfied */
852:       if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
853:           || ((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");
854:     }
855:     lf[i]      = want;
856:     startc    += lc[i];
857:     startf    += lf[i];
858:     remaining -= lf[i];
859:   }
860:   return(0);
861: }

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:     PetscMemcpy(lc,lf,m*sizeof(lf[0]));
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: }

908: PetscErrorCode  DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
909: {
911:   PetscInt       M,N,P,i;
912:   DM             da2;
913:   DM_DA          *dd = (DM_DA*)da->data,*dd2;


919:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
920:     M = dd->refine_x*dd->M;
921:   } else {
922:     M = 1 + dd->refine_x*(dd->M - 1);
923:   }
924:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
925:     if (dd->dim > 1) {
926:       N = dd->refine_y*dd->N;
927:     } else {
928:       N = 1;
929:     }
930:   } else {
931:     N = 1 + dd->refine_y*(dd->N - 1);
932:   }
933:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
934:     if (dd->dim > 2) {
935:       P = dd->refine_z*dd->P;
936:     } else {
937:       P = 1;
938:     }
939:   } else {
940:     P = 1 + dd->refine_z*(dd->P - 1);
941:   }
942:   DMDACreate(PetscObjectComm((PetscObject)da),&da2);
943:   DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
944:   DMDASetDim(da2,dd->dim);
945:   DMDASetSizes(da2,M,N,P);
946:   DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
947:   DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
948:   DMDASetDof(da2,dd->w);
949:   DMDASetStencilType(da2,dd->stencil_type);
950:   DMDASetStencilWidth(da2,dd->s);
951:   if (dd->dim == 3) {
952:     PetscInt *lx,*ly,*lz;
953:     PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
954:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
955:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
956:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
957:     DMDASetOwnershipRanges(da2,lx,ly,lz);
958:     PetscFree3(lx,ly,lz);
959:   } else if (dd->dim == 2) {
960:     PetscInt *lx,*ly;
961:     PetscMalloc2(dd->m,&lx,dd->n,&ly);
962:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
963:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
964:     DMDASetOwnershipRanges(da2,lx,ly,NULL);
965:     PetscFree2(lx,ly);
966:   } else if (dd->dim == 1) {
967:     PetscInt *lx;
968:     PetscMalloc1(dd->m,&lx);
969:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
970:     DMDASetOwnershipRanges(da2,lx,NULL,NULL);
971:     PetscFree(lx);
972:   }
973:   dd2 = (DM_DA*)da2->data;

975:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
976:   da2->ops->creatematrix = da->ops->creatematrix;
977:   /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
978:   da2->ops->getcoloring = da->ops->getcoloring;
979:   dd2->interptype       = dd->interptype;

981:   /* copy fill information if given */
982:   if (dd->dfill) {
983:     PetscMalloc1((dd->dfill[dd->w]+dd->w+1),&dd2->dfill);
984:     PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
985:   }
986:   if (dd->ofill) {
987:     PetscMalloc1((dd->ofill[dd->w]+dd->w+1),&dd2->ofill);
988:     PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
989:   }
990:   /* copy the refine information */
991:   dd2->coarsen_x = dd2->refine_x = dd->refine_x;
992:   dd2->coarsen_y = dd2->refine_y = dd->refine_y;
993:   dd2->coarsen_z = dd2->refine_z = dd->refine_z;

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

998:   dd2->lf = dd->lf;
999:   dd2->lj = dd->lj;

1001:   da2->leveldown = da->leveldown;
1002:   da2->levelup   = da->levelup + 1;

1004:   DMSetFromOptions(da2);
1005:   DMSetUp(da2);

1007:   /* interpolate coordinates if they are set on the coarse grid */
1008:   if (da->coordinates) {
1009:     DM  cdaf,cdac;
1010:     Vec coordsc,coordsf;
1011:     Mat II;

1013:     DMGetCoordinateDM(da,&cdac);
1014:     DMGetCoordinates(da,&coordsc);
1015:     DMGetCoordinateDM(da2,&cdaf);
1016:     /* force creation of the coordinate vector */
1017:     DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1018:     DMGetCoordinates(da2,&coordsf);
1019:     DMCreateInterpolation(cdac,cdaf,&II,NULL);
1020:     MatInterpolate(II,coordsc,coordsf);
1021:     MatDestroy(&II);
1022:   }

1024:   for (i=0; i<da->bs; i++) {
1025:     const char *fieldname;
1026:     DMDAGetFieldName(da,i,&fieldname);
1027:     DMDASetFieldName(da2,i,fieldname);
1028:   }

1030:   *daref = da2;
1031:   return(0);
1032: }


1037: PetscErrorCode  DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1038: {
1040:   PetscInt       M,N,P,i;
1041:   DM             da2;
1042:   DM_DA          *dd = (DM_DA*)da->data,*dd2;


1048:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1049:     M = dd->M / dd->coarsen_x;
1050:   } else {
1051:     M = 1 + (dd->M - 1) / dd->coarsen_x;
1052:   }
1053:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1054:     if (dd->dim > 1) {
1055:       N = dd->N / dd->coarsen_y;
1056:     } else {
1057:       N = 1;
1058:     }
1059:   } else {
1060:     N = 1 + (dd->N - 1) / dd->coarsen_y;
1061:   }
1062:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1063:     if (dd->dim > 2) {
1064:       P = dd->P / dd->coarsen_z;
1065:     } else {
1066:       P = 1;
1067:     }
1068:   } else {
1069:     P = 1 + (dd->P - 1) / dd->coarsen_z;
1070:   }
1071:   DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1072:   DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1073:   DMDASetDim(da2,dd->dim);
1074:   DMDASetSizes(da2,M,N,P);
1075:   DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1076:   DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1077:   DMDASetDof(da2,dd->w);
1078:   DMDASetStencilType(da2,dd->stencil_type);
1079:   DMDASetStencilWidth(da2,dd->s);
1080:   if (dd->dim == 3) {
1081:     PetscInt *lx,*ly,*lz;
1082:     PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1083:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1084:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1085:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1086:     DMDASetOwnershipRanges(da2,lx,ly,lz);
1087:     PetscFree3(lx,ly,lz);
1088:   } else if (dd->dim == 2) {
1089:     PetscInt *lx,*ly;
1090:     PetscMalloc2(dd->m,&lx,dd->n,&ly);
1091:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1092:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1093:     DMDASetOwnershipRanges(da2,lx,ly,NULL);
1094:     PetscFree2(lx,ly);
1095:   } else if (dd->dim == 1) {
1096:     PetscInt *lx;
1097:     PetscMalloc1(dd->m,&lx);
1098:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1099:     DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1100:     PetscFree(lx);
1101:   }
1102:   dd2 = (DM_DA*)da2->data;

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

1110:   /* copy fill information if given */
1111:   if (dd->dfill) {
1112:     PetscMalloc1((dd->dfill[dd->w]+dd->w+1),&dd2->dfill);
1113:     PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1114:   }
1115:   if (dd->ofill) {
1116:     PetscMalloc1((dd->ofill[dd->w]+dd->w+1),&dd2->ofill);
1117:     PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1118:   }
1119:   /* copy the refine information */
1120:   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1121:   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1122:   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;

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

1127:   dd2->lf = dd->lf;
1128:   dd2->lj = dd->lj;

1130:   da2->leveldown = da->leveldown + 1;
1131:   da2->levelup   = da->levelup;

1133:   DMSetFromOptions(da2);
1134:   DMSetUp(da2);

1136:   /* inject coordinates if they are set on the fine grid */
1137:   if (da->coordinates) {
1138:     DM         cdaf,cdac;
1139:     Vec        coordsc,coordsf;
1140:     VecScatter inject;

1142:     DMGetCoordinateDM(da,&cdaf);
1143:     DMGetCoordinates(da,&coordsf);
1144:     DMGetCoordinateDM(da2,&cdac);
1145:     /* force creation of the coordinate vector */
1146:     DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1147:     DMGetCoordinates(da2,&coordsc);

1149:     DMCreateInjection(cdac,cdaf,&inject);
1150:     VecScatterBegin(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1151:     VecScatterEnd(inject,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1152:     VecScatterDestroy(&inject);
1153:   }

1155:   for (i=0; i<da->bs; i++) {
1156:     const char *fieldname;
1157:     DMDAGetFieldName(da,i,&fieldname);
1158:     DMDASetFieldName(da2,i,fieldname);
1159:   }

1161:   *daref = da2;
1162:   return(0);
1163: }

1167: PetscErrorCode  DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1168: {
1170:   PetscInt       i,n,*refx,*refy,*refz;

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

1178:   /* Get refinement factors, defaults taken from the coarse DMDA */
1179:   PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1180:   for (i=0; i<nlevels; i++) {
1181:     DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1182:   }
1183:   n    = nlevels;
1184:   PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1185:   n    = nlevels;
1186:   PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1187:   n    = nlevels;
1188:   PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);

1190:   DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1191:   DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1192:   for (i=1; i<nlevels; i++) {
1193:     DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1194:     DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1195:   }
1196:   PetscFree3(refx,refy,refz);
1197:   return(0);
1198: }

1202: PetscErrorCode  DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1203: {
1205:   PetscInt       i;

1209:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1210:   if (nlevels == 0) return(0);
1212:   DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1213:   for (i=1; i<nlevels; i++) {
1214:     DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1215:   }
1216:   return(0);
1217: }