Actual source code: da.c

petsc-master 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/dmdaimpl.h>

  3: /*@
  4:   DMDASetSizes - Sets the number of grid points in the three dimensional directions

  6:   Logically Collective on da

  8:   Input Parameters:
  9: + da - the DMDA
 10: . M - the global X size
 11: . N - the global Y size
 12: - P - the global Z size

 14:   Level: intermediate

 16:   Developer Notes:
 17:   Since the dimension may not yet have been set the code cannot error check for non-positive Y and Z number of grid points

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

 30:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
 31:   if (M < 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in X direction must be positive");
 32:   if (N < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in Y direction must be positive");
 33:   if (P < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in Z direction must be positive");

 35:   dd->M = M;
 36:   dd->N = N;
 37:   dd->P = P;
 38:   return(0);
 39: }

 41: /*@
 42:   DMDASetNumProcs - Sets the number of processes in each dimension

 44:   Logically Collective on da

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

 52:   Level: intermediate

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

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

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

 88:   Not collective

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

 94:   Level: intermediate

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

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

114: /*@
115:   DMDASetDof - Sets the number of degrees of freedom per vertex

117:   Not collective

119:   Input Parameters:
120: + da  - The DMDA
121: - dof - Number of degrees of freedom

123:   Level: intermediate

125: .seealso: DMDAGetDof(), DMDACreate(), DMDestroy(), DMDA
126: @*/
127: PetscErrorCode  DMDASetDof(DM da, PetscInt dof)
128: {
129:   DM_DA *dd = (DM_DA*)da->data;

134:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
135:   dd->w  = dof;
136:   da->bs = dof;
137:   return(0);
138: }

140: /*@
141:   DMDAGetDof - Gets the number of degrees of freedom per vertex

143:   Not collective

145:   Input Parameter:
146: . da  - The DMDA

148:   Output Parameter:
149: . dof - Number of degrees of freedom

151:   Level: intermediate

153: .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
154: @*/
155: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
156: {
157:   DM_DA *dd = (DM_DA *) da->data;

162:   *dof = dd->w;
163:   return(0);
164: }

166: /*@
167:   DMDAGetOverlap - Gets the size of the per-processor overlap.

169:   Not collective

171:   Input Parameters:
172: . da  - The DMDA

174:   Output Parameters:
175: + x   - Overlap in the x direction
176: . y   - Overlap in the y direction
177: - z   - Overlap in the z direction

179:   Level: intermediate

181: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
182: @*/
183: PetscErrorCode  DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
184: {
185:   DM_DA *dd = (DM_DA*)da->data;

189:   if (x) *x = dd->xol;
190:   if (y) *y = dd->yol;
191:   if (z) *z = dd->zol;
192:   return(0);
193: }

195: /*@
196:   DMDASetOverlap - Sets the size of the per-processor overlap.

198:   Not collective

200:   Input Parameters:
201: + da  - The DMDA
202: . x   - Overlap in the x direction
203: . y   - Overlap in the y direction
204: - z   - Overlap in the z direction

206:   Level: intermediate

208: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
209: @*/
210: PetscErrorCode  DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
211: {
212:   DM_DA *dd = (DM_DA*)da->data;

219:   dd->xol = x;
220:   dd->yol = y;
221:   dd->zol = z;
222:   return(0);
223: }


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

229:   Not collective

231:   Input Parameters:
232: . da  - The DMDA

234:   Output Parameters:
235: + Nsub   - Number of local subdomains created upon decomposition

237:   Level: intermediate

239: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
240: @*/
241: PetscErrorCode  DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
242: {
243:   DM_DA *dd = (DM_DA*)da->data;

247:   if (Nsub) *Nsub = dd->Nsub;
248:   return(0);
249: }

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

254:   Not collective

256:   Input Parameters:
257: + da  - The DMDA
258: - Nsub - The number of local subdomains requested

260:   Level: intermediate

262: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
263: @*/
264: PetscErrorCode  DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
265: {
266:   DM_DA *dd = (DM_DA*)da->data;

271:   dd->Nsub = Nsub;
272:   return(0);
273: }

275: /*@
276:   DMDASetOffset - Sets the index offset of the DA.

278:   Collective on DA

280:   Input Parameter:
281: + da  - The DMDA
282: . xo  - The offset in the x direction
283: . yo  - The offset in the y direction
284: - zo  - The offset in the z direction

286:   Level: intermediate

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

292: .seealso: DMDAGetOffset(), DMDAVecGetArray()
293: @*/
294: PetscErrorCode  DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
295: {
297:   DM_DA          *dd = (DM_DA*)da->data;

307:   dd->xo = xo;
308:   dd->yo = yo;
309:   dd->zo = zo;
310:   dd->Mo = Mo;
311:   dd->No = No;
312:   dd->Po = Po;

314:   if (da->coordinateDM) {
315:     DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
316:   }
317:   return(0);
318: }

320: /*@
321:   DMDAGetOffset - Gets the index offset of the DA.

323:   Not collective

325:   Input Parameter:
326: . da  - The DMDA

328:   Output Parameters:
329: + xo  - The offset in the x direction
330: . yo  - The offset in the y direction
331: . zo  - The offset in the z direction
332: . Mo  - The global size in the x direction
333: . No  - The global size in the y direction
334: - Po  - The global size in the z direction

336:   Level: intermediate

338: .seealso: DMDASetOffset(), DMDAVecGetArray()
339: @*/
340: PetscErrorCode  DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
341: {
342:   DM_DA *dd = (DM_DA*)da->data;

346:   if (xo) *xo = dd->xo;
347:   if (yo) *yo = dd->yo;
348:   if (zo) *zo = dd->zo;
349:   if (Mo) *Mo = dd->Mo;
350:   if (No) *No = dd->No;
351:   if (Po) *Po = dd->Po;
352:   return(0);
353: }

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

358:   Not collective

360:   Input Parameter:
361: . da  - The DMDA

363:   Output Parameters:
364: + xs  - The start of the region in x
365: . ys  - The start of the region in y
366: . zs  - The start of the region in z
367: . xs  - The size of the region in x
368: . ys  - The size of the region in y
369: . zs  - The size of the region in z

371:   Level: intermediate

373: .seealso: DMDAGetOffset(), DMDAVecGetArray()
374: @*/
375: PetscErrorCode  DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
376: {
377:   DM_DA          *dd = (DM_DA*)da->data;

381:   if (xs) *xs = dd->nonxs;
382:   if (ys) *ys = dd->nonys;
383:   if (zs) *zs = dd->nonzs;
384:   if (xm) *xm = dd->nonxm;
385:   if (ym) *ym = dd->nonym;
386:   if (zm) *zm = dd->nonzm;
387:   return(0);
388: }


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

394:   Collective on DA

396:   Input Parameter:
397: + da  - The DMDA
398: . xs  - The start of the region in x
399: . ys  - The start of the region in y
400: . zs  - The start of the region in z
401: . xs  - The size of the region in x
402: . ys  - The size of the region in y
403: . zs  - The size of the region in z

405:   Level: intermediate

407: .seealso: DMDAGetOffset(), DMDAVecGetArray()
408: @*/
409: PetscErrorCode  DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
410: {
411:   DM_DA          *dd = (DM_DA*)da->data;

421:   dd->nonxs = xs;
422:   dd->nonys = ys;
423:   dd->nonzs = zs;
424:   dd->nonxm = xm;
425:   dd->nonym = ym;
426:   dd->nonzm = zm;

428:   return(0);
429: }

431: /*@
432:   DMDASetStencilType - Sets the type of the communication stencil

434:   Logically Collective on da

436:   Input Parameter:
437: + da    - The DMDA
438: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.

440:   Level: intermediate

442: .seealso: DMDACreate(), DMDestroy(), DMDA
443: @*/
444: PetscErrorCode  DMDASetStencilType(DM da, DMDAStencilType stype)
445: {
446:   DM_DA *dd = (DM_DA*)da->data;

451:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
452:   dd->stencil_type = stype;
453:   return(0);
454: }

456: /*@
457:   DMDAGetStencilType - Gets the type of the communication stencil

459:   Not collective

461:   Input Parameter:
462: . da    - The DMDA

464:   Output Parameter:
465: . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.

467:   Level: intermediate

469: .seealso: DMDACreate(), DMDestroy(), DMDA
470: @*/
471: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
472: {
473:   DM_DA *dd = (DM_DA*)da->data;

478:   *stype = dd->stencil_type;
479:   return(0);
480: }

482: /*@
483:   DMDASetStencilWidth - Sets the width of the communication stencil

485:   Logically Collective on da

487:   Input Parameter:
488: + da    - The DMDA
489: - width - The stencil width

491:   Level: intermediate

493: .seealso: DMDACreate(), DMDestroy(), DMDA
494: @*/
495: PetscErrorCode  DMDASetStencilWidth(DM da, PetscInt width)
496: {
497:   DM_DA *dd = (DM_DA*)da->data;

502:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
503:   dd->s = width;
504:   return(0);
505: }

507: /*@
508:   DMDAGetStencilWidth - Gets the width of the communication stencil

510:   Not collective

512:   Input Parameter:
513: . da    - The DMDA

515:   Output Parameter:
516: . width - The stencil width

518:   Level: intermediate

520: .seealso: DMDACreate(), DMDestroy(), DMDA
521: @*/
522: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
523: {
524:   DM_DA *dd = (DM_DA *) da->data;

529:   *width = dd->s;
530:   return(0);
531: }

533: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
534: {
535:   PetscInt i,sum;

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

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

547:   Logically Collective on da

549:   Input Parameter:
550: + da - The DMDA
551: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
552: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
553: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.

555:   Level: intermediate

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

559: .seealso: DMDACreate(), DMDestroy(), DMDA
560: @*/
561: PetscErrorCode  DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
562: {
564:   DM_DA          *dd = (DM_DA*)da->data;

568:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
569:   if (lx) {
570:     if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
571:     DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
572:     if (!dd->lx) {
573:       PetscMalloc1(dd->m, &dd->lx);
574:     }
575:     PetscArraycpy(dd->lx, lx, dd->m);
576:   }
577:   if (ly) {
578:     if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
579:     DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
580:     if (!dd->ly) {
581:       PetscMalloc1(dd->n, &dd->ly);
582:     }
583:     PetscArraycpy(dd->ly, ly, dd->n);
584:   }
585:   if (lz) {
586:     if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
587:     DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
588:     if (!dd->lz) {
589:       PetscMalloc1(dd->p, &dd->lz);
590:     }
591:     PetscArraycpy(dd->lz, lz, dd->p);
592:   }
593:   return(0);
594: }

596: /*@
597:        DMDASetInterpolationType - Sets the type of interpolation that will be
598:           returned by DMCreateInterpolation()

600:    Logically Collective on da

602:    Input Parameter:
603: +  da - initial distributed array
604: .  ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms

606:    Level: intermediate

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

611: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
612: @*/
613: PetscErrorCode  DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
614: {
615:   DM_DA *dd = (DM_DA*)da->data;

620:   dd->interptype = ctype;
621:   return(0);
622: }

624: /*@
625:        DMDAGetInterpolationType - Gets the type of interpolation that will be
626:           used by DMCreateInterpolation()

628:    Not Collective

630:    Input Parameter:
631: .  da - distributed array

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

636:    Level: intermediate

638: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
639: @*/
640: PetscErrorCode  DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
641: {
642:   DM_DA *dd = (DM_DA*)da->data;

647:   *ctype = dd->interptype;
648:   return(0);
649: }

651: /*@C
652:       DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
653:         processes neighbors.

655:     Not Collective

657:    Input Parameter:
658: .     da - the DMDA object

660:    Output Parameters:
661: .     ranks - the neighbors ranks, stored with the x index increasing most rapidly.
662:               this process itself is in the list

664:    Notes:
665:     In 2d the array is of length 9, in 3d of length 27
666:           Not supported in 1d
667:           Do not free the array, it is freed when the DMDA is destroyed.

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

672:    Level: intermediate

674: @*/
675: PetscErrorCode  DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
676: {
677:   DM_DA *dd = (DM_DA*)da->data;

681:   *ranks = dd->neighbors;
682:   return(0);
683: }

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

688:     Not Collective

690:    Input Parameter:
691: .     da - the DMDA object

693:    Output Parameter:
694: +     lx - ownership along x direction (optional)
695: .     ly - ownership along y direction (optional)
696: -     lz - ownership along z direction (optional)

698:    Level: intermediate

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

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

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

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

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

718:   if (lx) *lx = dd->lx;
719:   if (ly) *ly = dd->ly;
720:   if (lz) *lz = dd->lz;
721:   return(0);
722: }

724: /*@
725:      DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined

727:     Logically Collective on da

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

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

740:   Level: intermediate

742:     Notes:
743:     Pass PETSC_IGNORE to leave a value unchanged

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


757:   if (refine_x > 0) dd->refine_x = refine_x;
758:   if (refine_y > 0) dd->refine_y = refine_y;
759:   if (refine_z > 0) dd->refine_z = refine_z;
760:   return(0);
761: }

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

766:     Not Collective

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

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

776:   Level: intermediate

778:     Notes:
779:     Pass NULL for values you do not need

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

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

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

798:     Logically Collective on da

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:
807:     See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
808:        the diagonal and off-diagonal blocks of the matrix

810:    Not supported from Fortran

812: .seealso: DMCreateMatrix(), DMDASetBlockFills()
813: @*/
814: PetscErrorCode  DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
815: {
818:   da->ops->creatematrix = f;
819:   return(0);
820: }

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

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

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

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

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

876:   if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
877:   if (ratio == 1) {
878:     PetscArraycpy(lc,lf,m);
879:     return(0);
880:   }
881:   for (i=0,totalf=0; i<m; i++) totalf += lf[i];
882:   remaining = (!periodic) + (totalf - (!periodic)) / ratio;
883:   for (i=0,startc=0,startf=0; i<m; i++) {
884:     PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
885:     if (i == m-1) lc[i] = want;
886:     else {
887:       const PetscInt nextf = startf+lf[i];
888:       /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
889:        * node is within one stencil width. */
890:       while (nextf/ratio < startc+want-stencil_width) want--;
891:       /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
892:        * fine node is within one stencil width. */
893:       while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
894:       if (want < 0 || want > remaining
895:           || (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");
896:     }
897:     lc[i]      = want;
898:     startc    += lc[i];
899:     startf    += lf[i];
900:     remaining -= lc[i];
901:   }
902:   return(0);
903: }

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


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

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

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

993:   if (dd->refine_z_hier) {
994:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) {
995:       dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
996:     }
997:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) {
998:       dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
999:     }
1000:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1001:     PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1002:     PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);
1003:   }
1004:   if (dd->refine_y_hier) {
1005:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) {
1006:       dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
1007:     }
1008:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) {
1009:       dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
1010:     }
1011:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1012:     PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1013:     PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);
1014:   }
1015:   if (dd->refine_x_hier) {
1016:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) {
1017:       dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
1018:     }
1019:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) {
1020:       dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
1021:     }
1022:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1023:     PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1024:     PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);
1025:   }


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

1031:   dd2->lf = dd->lf;
1032:   dd2->lj = dd->lj;

1034:   da2->leveldown = da->leveldown;
1035:   da2->levelup   = da->levelup + 1;

1037:   DMSetUp(da2);

1039:   /* interpolate coordinates if they are set on the coarse grid */
1040:   if (da->coordinates) {
1041:     DM  cdaf,cdac;
1042:     Vec coordsc,coordsf;
1043:     Mat II;

1045:     DMGetCoordinateDM(da,&cdac);
1046:     DMGetCoordinates(da,&coordsc);
1047:     DMGetCoordinateDM(da2,&cdaf);
1048:     /* force creation of the coordinate vector */
1049:     DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1050:     DMGetCoordinates(da2,&coordsf);
1051:     DMCreateInterpolation(cdac,cdaf,&II,NULL);
1052:     MatInterpolate(II,coordsc,coordsf);
1053:     MatDestroy(&II);
1054:   }

1056:   for (i=0; i<da->bs; i++) {
1057:     const char *fieldname;
1058:     DMDAGetFieldName(da,i,&fieldname);
1059:     DMDASetFieldName(da2,i,fieldname);
1060:   }

1062:   *daref = da2;
1063:   return(0);
1064: }


1067: PetscErrorCode  DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1068: {
1070:   PetscInt       M,N,P,i,dim;
1071:   DM             da2;
1072:   DM_DA          *dd = (DM_DA*)da->data,*dd2;


1078:   DMGetDimension(da, &dim);
1079:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1080:     M = dd->M / dd->coarsen_x;
1081:   } else {
1082:     M = 1 + (dd->M - 1) / dd->coarsen_x;
1083:   }
1084:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1085:     if (dim > 1) {
1086:       N = dd->N / dd->coarsen_y;
1087:     } else {
1088:       N = 1;
1089:     }
1090:   } else {
1091:     N = 1 + (dd->N - 1) / dd->coarsen_y;
1092:   }
1093:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1094:     if (dim > 2) {
1095:       P = dd->P / dd->coarsen_z;
1096:     } else {
1097:       P = 1;
1098:     }
1099:   } else {
1100:     P = 1 + (dd->P - 1) / dd->coarsen_z;
1101:   }
1102:   DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1103:   DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1104:   DMSetDimension(da2,dim);
1105:   DMDASetSizes(da2,M,N,P);
1106:   DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1107:   DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1108:   DMDASetDof(da2,dd->w);
1109:   DMDASetStencilType(da2,dd->stencil_type);
1110:   DMDASetStencilWidth(da2,dd->s);
1111:   if (dim == 3) {
1112:     PetscInt *lx,*ly,*lz;
1113:     PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1114:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1115:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1116:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1117:     DMDASetOwnershipRanges(da2,lx,ly,lz);
1118:     PetscFree3(lx,ly,lz);
1119:   } else if (dim == 2) {
1120:     PetscInt *lx,*ly;
1121:     PetscMalloc2(dd->m,&lx,dd->n,&ly);
1122:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1123:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1124:     DMDASetOwnershipRanges(da2,lx,ly,NULL);
1125:     PetscFree2(lx,ly);
1126:   } else if (dim == 1) {
1127:     PetscInt *lx;
1128:     PetscMalloc1(dd->m,&lx);
1129:     DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1130:     DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1131:     PetscFree(lx);
1132:   }
1133:   dd2 = (DM_DA*)da2->data;

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

1141:   /* copy fill information if given */
1142:   if (dd->dfill) {
1143:     PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1144:     PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);
1145:   }
1146:   if (dd->ofill) {
1147:     PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1148:     PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);
1149:   }
1150:   /* copy the refine information */
1151:   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1152:   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1153:   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;

1155:   if (dd->refine_z_hier) {
1156:     if (da->levelup - da->leveldown -1 > -1 && da->levelup - da->leveldown - 1< dd->refine_z_hier_n) {
1157:       dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown - 1];
1158:     }
1159:     if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_z_hier_n) {
1160:       dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown - 2];
1161:     }
1162:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1163:     PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1164:     PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);
1165:   }
1166:   if (dd->refine_y_hier) {
1167:     if (da->levelup - da->leveldown - 1 > -1 && da->levelup - da->leveldown - 1< dd->refine_y_hier_n) {
1168:       dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown - 1];
1169:     }
1170:     if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_y_hier_n) {
1171:       dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown - 2];
1172:     }
1173:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1174:     PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1175:     PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);
1176:   }
1177:   if (dd->refine_x_hier) {
1178:     if (da->levelup - da->leveldown - 1 > -1 && da->levelup - da->leveldown - 1 < dd->refine_x_hier_n) {
1179:       dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown - 1];
1180:     }
1181:     if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_x_hier_n) {
1182:       dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown - 2];
1183:     }
1184:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1185:     PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1186:     PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);
1187:   }

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

1192:   dd2->lf = dd->lf;
1193:   dd2->lj = dd->lj;

1195:   da2->leveldown = da->leveldown + 1;
1196:   da2->levelup   = da->levelup;

1198:   DMSetUp(da2);

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

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

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

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

1227:   *daref = da2;
1228:   return(0);
1229: }

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: }

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

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

1281: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm,PetscInt n,PetscReal *nodes)
1282: {
1284:   PetscInt       i,j,xs,xn,q;
1285:   PetscScalar    *xx;
1286:   PetscReal      h;
1287:   Vec            x;
1288:   DM_DA          *da = (DM_DA*)dm->data;

1291:   if (da->bx != DM_BOUNDARY_PERIODIC) {
1292:     DMDAGetInfo(dm,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1293:     q    = (q-1)/(n-1);  /* number of spectral elements */
1294:     h    = 2.0/q;
1295:     DMDAGetCorners(dm,&xs,NULL,NULL,&xn,NULL,NULL);
1296:     xs   = xs/(n-1);
1297:     xn   = xn/(n-1);
1298:     DMDASetUniformCoordinates(dm,-1.,1.,0.,0.,0.,0.);
1299:     DMGetCoordinates(dm,&x);
1300:     DMDAVecGetArray(dm,x,&xx);

1302:     /* loop over local spectral elements */
1303:     for (j=xs; j<xs+xn; j++) {
1304:       /*
1305:        Except for the first process, each process starts on the second GLL point of the first element on that process
1306:        */
1307:       for (i= (j == xs && xs > 0)? 1 : 0; i<n; i++) {
1308:         xx[j*(n-1) + i] = -1.0 + h*j + h*(nodes[i]+1.0)/2.;
1309:       }
1310:     }
1311:     DMDAVecRestoreArray(dm,x,&xx);
1312:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for periodic");
1313:   return(0);
1314: }

1316: /*@

1318:      DMDASetGLLCoordinates - Sets the global coordinates from -1 to 1 to the GLL points of as many GLL elements that fit the number of grid points

1320:    Collective on da

1322:    Input Parameters:
1323: +   da - the DMDA object
1324: -   n - the number of GLL nodes
1325: -   nodes - the GLL nodes

1327:    Notes:
1328:     the parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1329:           on each process much be divisible by the number of GLL elements needed per process. This depends on whether the DM is
1330:           periodic or not.

1332:    Level: advanced

1334: .seealso:   DMDACreate(), PetscDTGaussLobattoLegendreQuadrature(), DMGetCoordinates()
1335: @*/
1336: PetscErrorCode DMDASetGLLCoordinates(DM da,PetscInt n,PetscReal *nodes)
1337: {

1341:   if (da->dim == 1) {
1342:     DMDASetGLLCoordinates_1d(da,n,nodes);
1343:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for 2 or 3d");
1344:   return(0);
1345: }

1347: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM da1,DM dm2,PetscBool *compatible,PetscBool *set)
1348: {
1350:   DM_DA          *dd1 = (DM_DA*)da1->data,*dd2;
1351:   DM             da2;
1352:   DMType         dmtype2;
1353:   PetscBool      isda,compatibleLocal;
1354:   PetscInt       i;

1357:   if (!da1->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da1),PETSC_ERR_ARG_WRONGSTATE,"DMSetUp() must be called on first DM before DMGetCompatibility()");
1358:   DMGetType(dm2,&dmtype2);
1359:   PetscStrcmp(dmtype2,DMDA,&isda);
1360:   if (isda) {
1361:     da2 = dm2;
1362:     dd2 = (DM_DA*)da2->data;
1363:     if (!da2->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da2),PETSC_ERR_ARG_WRONGSTATE,"DMSetUp() must be called on second DM before DMGetCompatibility()");
1364:     compatibleLocal = (PetscBool)(da1->dim == da2->dim);
1365:     if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
1366:     /*                                                                           Global size              ranks               Boundary type */
1367:     if (compatibleLocal)                 compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
1368:     if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
1369:     if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
1370:     if (compatibleLocal) {
1371:       for (i=0; i<dd1->m; ++i) {
1372:         compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i]));           /* Local size     */
1373:       }
1374:     }
1375:     if (compatibleLocal && da1->dim > 1) {
1376:       for (i=0; i<dd1->n; ++i) {
1377:         compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
1378:       }
1379:     }
1380:     if (compatibleLocal && da1->dim > 2) {
1381:       for (i=0; i<dd1->p; ++i) {
1382:         compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
1383:       }
1384:     }
1385:     *compatible = compatibleLocal;
1386:     *set = PETSC_TRUE;
1387:   } else {
1388:     /* Decline to determine compatibility with other DM types */
1389:     *set = PETSC_FALSE;
1390:   }
1391:   return(0);
1392: }