Actual source code: da.c

petsc-3.9.0 2018-04-07
Report Typos and Errors
  1:  #include <petsc/private/dmdaimpl.h>

  3: /*@
  4:   DMDASetSizes - Sets the global sizes

  6:   Logically Collective on DMDA

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

 14:   Level: intermediate

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

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

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

 35: /*@
 36:   DMDASetNumProcs - Sets the number of processes in each dimension

 38:   Logically Collective on DMDA

 40:   Input Parameters:
 41: + da - the DMDA
 42: . m - the number of X procs (or PETSC_DECIDE)
 43: . n - the number of Y procs (or PETSC_DECIDE)
 44: - p - the number of Z procs (or PETSC_DECIDE)

 46:   Level: intermediate

 48: .seealso: DMDASetSizes(), DMDAGetSize(), PetscSplitOwnership()
 49: @*/
 50: PetscErrorCode  DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
 51: {
 52:   DM_DA          *dd = (DM_DA*)da->data;

 60:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
 61:   dd->m = m;
 62:   dd->n = n;
 63:   dd->p = p;
 64:   if (da->dim == 2) {
 65:     PetscMPIInt size;
 66:     MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
 67:     if ((dd->m > 0) && (dd->n < 0)) {
 68:       dd->n = size/dd->m;
 69:       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);
 70:     }
 71:     if ((dd->n > 0) && (dd->m < 0)) {
 72:       dd->m = size/dd->n;
 73:       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);
 74:     }
 75:   }
 76:   return(0);
 77: }

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

 82:   Not collective

 84:   Input Parameter:
 85: + da    - The DMDA
 86: - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC

 88:   Level: intermediate

 90: .keywords:  distributed array, periodicity
 91: .seealso: DMDACreate(), DMDestroy(), DMDA, DMBoundaryType
 92: @*/
 93: PetscErrorCode  DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)
 94: {
 95:   DM_DA *dd = (DM_DA*)da->data;

102:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
103:   dd->bx = bx;
104:   dd->by = by;
105:   dd->bz = bz;
106:   return(0);
107: }

109: /*@
110:   DMDASetDof - Sets the number of degrees of freedom per vertex

112:   Not collective

114:   Input Parameters:
115: + da  - The DMDA
116: - dof - Number of degrees of freedom

118:   Level: intermediate

120: .keywords:  distributed array, degrees of freedom
121: .seealso: DMDAGetDof(), DMDACreate(), DMDestroy(), DMDA
122: @*/
123: PetscErrorCode  DMDASetDof(DM da, PetscInt dof)
124: {
125:   DM_DA *dd = (DM_DA*)da->data;

130:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
131:   dd->w  = dof;
132:   da->bs = dof;
133:   return(0);
134: }

136: /*@
137:   DMDAGetDof - Gets the number of degrees of freedom per vertex

139:   Not collective

141:   Input Parameter:
142: . da  - The DMDA

144:   Output Parameter:
145: . dof - Number of degrees of freedom

147:   Level: intermediate

149: .keywords:  distributed array, degrees of freedom
150: .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
151: @*/
152: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
153: {
154:   DM_DA *dd = (DM_DA *) da->data;

159:   *dof = dd->w;
160:   return(0);
161: }

163: /*@
164:   DMDAGetOverlap - Gets the size of the per-processor overlap.

166:   Not collective

168:   Input Parameters:
169: . da  - The DMDA

171:   Output Parameters:
172: + x   - Overlap in the x direction
173: . y   - Overlap in the y direction
174: - z   - Overlap in the z direction

176:   Level: intermediate

178: .keywords:  distributed array, overlap, domain decomposition
179: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
180: @*/
181: PetscErrorCode  DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
182: {
183:   DM_DA *dd = (DM_DA*)da->data;

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

193: /*@
194:   DMDASetOverlap - Sets the size of the per-processor overlap.

196:   Not collective

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

204:   Level: intermediate

206: .keywords:  distributed array, overlap, domain decomposition
207: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
208: @*/
209: PetscErrorCode  DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
210: {
211:   DM_DA *dd = (DM_DA*)da->data;

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


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

228:   Not collective

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

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

236:   Level: intermediate

238: .keywords:  distributed array, domain decomposition
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: .keywords:  distributed array, domain decomposition
263: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
264: @*/
265: PetscErrorCode  DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
266: {
267:   DM_DA *dd = (DM_DA*)da->data;

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

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

279:   Collective on DA

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

287:   Level: intermediate

289:   Notes: 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: .keywords:  distributed array, degrees of freedom
293: .seealso: DMDAGetOffset(), DMDAVecGetArray()
294: @*/
295: PetscErrorCode  DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
296: {
298:   DM_DA          *dd = (DM_DA*)da->data;

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

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

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

324:   Not collective

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

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

337:   Level: intermediate

339: .keywords:  distributed array, degrees of freedom
340: .seealso: DMDASetOffset(), DMDAVecGetArray()
341: @*/
342: PetscErrorCode  DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
343: {
344:   DM_DA *dd = (DM_DA*)da->data;

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

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

360:   Not collective

362:   Input Parameter:
363: . da  - The DMDA

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

373:   Level: intermediate

375: .keywords:  distributed array, degrees of freedom
376: .seealso: DMDAGetOffset(), DMDAVecGetArray()
377: @*/
378: PetscErrorCode  DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
379: {
380:   DM_DA          *dd = (DM_DA*)da->data;

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


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

397:   Collective on DA

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

408:   Level: intermediate

410: .keywords:  distributed array, degrees of freedom
411: .seealso: DMDAGetOffset(), DMDAVecGetArray()
412: @*/
413: PetscErrorCode  DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
414: {
415:   DM_DA          *dd = (DM_DA*)da->data;

425:   dd->nonxs = xs;
426:   dd->nonys = ys;
427:   dd->nonzs = zs;
428:   dd->nonxm = xm;
429:   dd->nonym = ym;
430:   dd->nonzm = zm;

432:   return(0);
433: }

435: /*@
436:   DMDASetStencilType - Sets the type of the communication stencil

438:   Logically Collective on DMDA

440:   Input Parameter:
441: + da    - The DMDA
442: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.

444:   Level: intermediate

446: .keywords:  distributed array, stencil
447: .seealso: DMDACreate(), DMDestroy(), DMDA
448: @*/
449: PetscErrorCode  DMDASetStencilType(DM da, DMDAStencilType stype)
450: {
451:   DM_DA *dd = (DM_DA*)da->data;

456:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
457:   dd->stencil_type = stype;
458:   return(0);
459: }

461: /*@
462:   DMDAGetStencilType - Gets the type of the communication stencil

464:   Not collective

466:   Input Parameter:
467: . da    - The DMDA

469:   Output Parameter:
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 DMDAGetStencilType(DM da, DMDAStencilType *stype)
478: {
479:   DM_DA *dd = (DM_DA*)da->data;

484:   *stype = dd->stencil_type;
485:   return(0);
486: }

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

491:   Logically Collective on DMDA

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

497:   Level: intermediate

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

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

514: /*@
515:   DMDAGetStencilWidth - Gets the width of the communication stencil

517:   Not collective

519:   Input Parameter:
520: . da    - The DMDA

522:   Output Parameter:
523: . width - The stencil width

525:   Level: intermediate

527: .keywords:  distributed array, stencil
528: .seealso: DMDACreate(), DMDestroy(), DMDA
529: @*/
530: PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
531: {
532:   DM_DA *dd = (DM_DA *) da->data;

537:   *width = dd->s;
538:   return(0);
539: }

541: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
542: {
543:   PetscInt i,sum;

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

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

555:   Logically Collective on DMDA

557:   Input Parameter:
558: + da - The DMDA
559: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
560: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
561: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.

563:   Level: intermediate

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

567: .keywords:  distributed array
568: .seealso: DMDACreate(), DMDestroy(), DMDA
569: @*/
570: PetscErrorCode  DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
571: {
573:   DM_DA          *dd = (DM_DA*)da->data;

577:   if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
578:   if (lx) {
579:     if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
580:     DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
581:     if (!dd->lx) {
582:       PetscMalloc1(dd->m, &dd->lx);
583:     }
584:     PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
585:   }
586:   if (ly) {
587:     if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
588:     DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
589:     if (!dd->ly) {
590:       PetscMalloc1(dd->n, &dd->ly);
591:     }
592:     PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
593:   }
594:   if (lz) {
595:     if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
596:     DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
597:     if (!dd->lz) {
598:       PetscMalloc1(dd->p, &dd->lz);
599:     }
600:     PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
601:   }
602:   return(0);
603: }

605: /*@
606:        DMDASetInterpolationType - Sets the type of interpolation that will be
607:           returned by DMCreateInterpolation()

609:    Logically Collective on DMDA

611:    Input Parameter:
612: +  da - initial distributed array
613: .  ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms

615:    Level: intermediate

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

619: .keywords:  distributed array, interpolation

621: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
622: @*/
623: PetscErrorCode  DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
624: {
625:   DM_DA *dd = (DM_DA*)da->data;

630:   dd->interptype = ctype;
631:   return(0);
632: }

634: /*@
635:        DMDAGetInterpolationType - Gets the type of interpolation that will be
636:           used by DMCreateInterpolation()

638:    Not Collective

640:    Input Parameter:
641: .  da - distributed array

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

646:    Level: intermediate

648: .keywords:  distributed array, interpolation

650: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
651: @*/
652: PetscErrorCode  DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
653: {
654:   DM_DA *dd = (DM_DA*)da->data;

659:   *ctype = dd->interptype;
660:   return(0);
661: }

663: /*@C
664:       DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
665:         processes neighbors.

667:     Not Collective

669:    Input Parameter:
670: .     da - the DMDA object

672:    Output Parameters:
673: .     ranks - the neighbors ranks, stored with the x index increasing most rapidly.
674:               this process itself is in the list

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

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

682:    Level: intermediate

684: @*/
685: PetscErrorCode  DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
686: {
687:   DM_DA *dd = (DM_DA*)da->data;

691:   *ranks = dd->neighbors;
692:   return(0);
693: }

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

698:     Not Collective

700:    Input Parameter:
701: .     da - the DMDA object

703:    Output Parameter:
704: +     lx - ownership along x direction (optional)
705: .     ly - ownership along y direction (optional)
706: -     lz - ownership along z direction (optional)

708:    Level: intermediate

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

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

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

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

720:      Not available from Fortran

722: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
723: @*/
724: PetscErrorCode  DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
725: {
726:   DM_DA *dd = (DM_DA*)da->data;

730:   if (lx) *lx = dd->lx;
731:   if (ly) *ly = dd->ly;
732:   if (lz) *lz = dd->lz;
733:   return(0);
734: }

736: /*@
737:      DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined

739:     Logically Collective on DMDA

741:   Input Parameters:
742: +    da - the DMDA object
743: .    refine_x - ratio of fine grid to coarse in x direction (2 by default)
744: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
745: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

747:   Options Database:
748: +  -da_refine_x - refinement ratio in x direction
749: .  -da_refine_y - refinement ratio in y direction
750: -  -da_refine_z - refinement ratio in z direction

752:   Level: intermediate

754:     Notes: Pass PETSC_IGNORE to leave a value unchanged

756: .seealso: DMRefine(), DMDAGetRefinementFactor()
757: @*/
758: PetscErrorCode  DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
759: {
760:   DM_DA *dd = (DM_DA*)da->data;


768:   if (refine_x > 0) dd->refine_x = refine_x;
769:   if (refine_y > 0) dd->refine_y = refine_y;
770:   if (refine_z > 0) dd->refine_z = refine_z;
771:   return(0);
772: }

774: /*@C
775:      DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined

777:     Not Collective

779:   Input Parameter:
780: .    da - the DMDA object

782:   Output Parameters:
783: +    refine_x - ratio of fine grid to coarse in x direction (2 by default)
784: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
785: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

787:   Level: intermediate

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

791: .seealso: DMRefine(), DMDASetRefinementFactor()
792: @*/
793: PetscErrorCode  DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
794: {
795:   DM_DA *dd = (DM_DA*)da->data;

799:   if (refine_x) *refine_x = dd->refine_x;
800:   if (refine_y) *refine_y = dd->refine_y;
801:   if (refine_z) *refine_z = dd->refine_z;
802:   return(0);
803: }

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

808:     Logically Collective on DMDA

810:   Input Parameters:
811: +    da - the DMDA object
812: -    f - the function that allocates the matrix for that specific DMDA

814:   Level: developer

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

819:    Not supported from Fortran

821: .seealso: DMCreateMatrix(), DMDASetBlockFills()
822: @*/
823: PetscErrorCode  DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
824: {
827:   da->ops->creatematrix = f;
828:   return(0);
829: }

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

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

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

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

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

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

914: PetscErrorCode  DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
915: {
917:   PetscInt       M,N,P,i,dim;
918:   DM             da2;
919:   DM_DA          *dd = (DM_DA*)da->data,*dd2;


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

982:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
983:   da2->ops->creatematrix = da->ops->creatematrix;
984:   /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
985:   da2->ops->getcoloring = da->ops->getcoloring;
986:   dd2->interptype       = dd->interptype;

988:   /* copy fill information if given */
989:   if (dd->dfill) {
990:     PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
991:     PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
992:   }
993:   if (dd->ofill) {
994:     PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
995:     PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
996:   }
997:   /* copy the refine information */
998:   dd2->coarsen_x = dd2->refine_x = dd->refine_x;
999:   dd2->coarsen_y = dd2->refine_y = dd->refine_y;
1000:   dd2->coarsen_z = dd2->refine_z = dd->refine_z;

1002:   if (dd->refine_z_hier) {
1003:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) {
1004:       dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
1005:     }
1006:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) {
1007:       dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
1008:     }
1009:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1010:     PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1011:     PetscMemcpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n*sizeof(PetscInt));
1012:   }
1013:   if (dd->refine_y_hier) {
1014:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) {
1015:       dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
1016:     }
1017:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) {
1018:       dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
1019:     }
1020:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1021:     PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1022:     PetscMemcpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n*sizeof(PetscInt));
1023:   }
1024:   if (dd->refine_x_hier) {
1025:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) {
1026:       dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
1027:     }
1028:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) {
1029:       dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
1030:     }
1031:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1032:     PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1033:     PetscMemcpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n*sizeof(PetscInt));
1034:   }
1035: 

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

1040:   dd2->lf = dd->lf;
1041:   dd2->lj = dd->lj;

1043:   da2->leveldown = da->leveldown;
1044:   da2->levelup   = da->levelup + 1;

1046:   DMSetUp(da2);

1048:   /* interpolate coordinates if they are set on the coarse grid */
1049:   if (da->coordinates) {
1050:     DM  cdaf,cdac;
1051:     Vec coordsc,coordsf;
1052:     Mat II;

1054:     DMGetCoordinateDM(da,&cdac);
1055:     DMGetCoordinates(da,&coordsc);
1056:     DMGetCoordinateDM(da2,&cdaf);
1057:     /* force creation of the coordinate vector */
1058:     DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1059:     DMGetCoordinates(da2,&coordsf);
1060:     DMCreateInterpolation(cdac,cdaf,&II,NULL);
1061:     MatInterpolate(II,coordsc,coordsf);
1062:     MatDestroy(&II);
1063:   }

1065:   for (i=0; i<da->bs; i++) {
1066:     const char *fieldname;
1067:     DMDAGetFieldName(da,i,&fieldname);
1068:     DMDASetFieldName(da2,i,fieldname);
1069:   }

1071:   *daref = da2;
1072:   return(0);
1073: }


1076: PetscErrorCode  DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1077: {
1079:   PetscInt       M,N,P,i,dim;
1080:   DM             da2;
1081:   DM_DA          *dd = (DM_DA*)da->data,*dd2;


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

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

1150:   /* copy fill information if given */
1151:   if (dd->dfill) {
1152:     PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1153:     PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1154:   }
1155:   if (dd->ofill) {
1156:     PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1157:     PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1158:   }
1159:   /* copy the refine information */
1160:   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1161:   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1162:   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;

1164:   if (dd->refine_z_hier) {
1165:     if (da->levelup - da->leveldown -1 > -1 && da->levelup - da->leveldown - 1< dd->refine_z_hier_n) {
1166:       dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown - 1];
1167:     }
1168:     if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_z_hier_n) {
1169:       dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown - 2];
1170:     }
1171:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1172:     PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1173:     PetscMemcpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n*sizeof(PetscInt));
1174:   }
1175:   if (dd->refine_y_hier) {
1176:     if (da->levelup - da->leveldown - 1 > -1 && da->levelup - da->leveldown - 1< dd->refine_y_hier_n) {
1177:       dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown - 1];
1178:     }
1179:     if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_y_hier_n) {
1180:       dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown - 2];
1181:     }
1182:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1183:     PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1184:     PetscMemcpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n*sizeof(PetscInt));
1185:   }
1186:   if (dd->refine_x_hier) {
1187:     if (da->levelup - da->leveldown - 1 > -1 && da->levelup - da->leveldown - 1 < dd->refine_x_hier_n) {
1188:       dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown - 1];
1189:     }
1190:     if (da->levelup - da->leveldown - 2 > -1 && da->levelup - da->leveldown - 2 < dd->refine_x_hier_n) {
1191:       dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown - 2];
1192:     }
1193:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1194:     PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1195:     PetscMemcpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n*sizeof(PetscInt));
1196:   }

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

1201:   dd2->lf = dd->lf;
1202:   dd2->lj = dd->lj;

1204:   da2->leveldown = da->leveldown + 1;
1205:   da2->levelup   = da->levelup;

1207:   DMSetUp(da2);

1209:   /* inject coordinates if they are set on the fine grid */
1210:   if (da->coordinates) {
1211:     DM         cdaf,cdac;
1212:     Vec        coordsc,coordsf;
1213:     Mat        inject;
1214:     VecScatter vscat;

1216:     DMGetCoordinateDM(da,&cdaf);
1217:     DMGetCoordinates(da,&coordsf);
1218:     DMGetCoordinateDM(da2,&cdac);
1219:     /* force creation of the coordinate vector */
1220:     DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1221:     DMGetCoordinates(da2,&coordsc);

1223:     DMCreateInjection(cdac,cdaf,&inject);
1224:     MatScatterGetVecScatter(inject,&vscat);
1225:     VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1226:     VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1227:     MatDestroy(&inject);
1228:   }

1230:   for (i=0; i<da->bs; i++) {
1231:     const char *fieldname;
1232:     DMDAGetFieldName(da,i,&fieldname);
1233:     DMDASetFieldName(da2,i,fieldname);
1234:   }

1236:   *daref = da2;
1237:   return(0);
1238: }

1240: PetscErrorCode  DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1241: {
1243:   PetscInt       i,n,*refx,*refy,*refz;

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

1251:   /* Get refinement factors, defaults taken from the coarse DMDA */
1252:   PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1253:   for (i=0; i<nlevels; i++) {
1254:     DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1255:   }
1256:   n    = nlevels;
1257:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1258:   n    = nlevels;
1259:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1260:   n    = nlevels;
1261:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);

1263:   DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1264:   DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1265:   for (i=1; i<nlevels; i++) {
1266:     DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1267:     DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1268:   }
1269:   PetscFree3(refx,refy,refz);
1270:   return(0);
1271: }

1273: PetscErrorCode  DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1274: {
1276:   PetscInt       i;

1280:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1281:   if (nlevels == 0) return(0);
1283:   DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1284:   for (i=1; i<nlevels; i++) {
1285:     DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1286:   }
1287:   return(0);
1288: }

1290:  #include <petscgll.h>

1292: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm,PetscGLL *gll)
1293: {
1295:   PetscInt       i,j,n = gll->n,xs,xn,q;
1296:   PetscScalar    *xx;
1297:   PetscReal      h;
1298:   Vec            x;
1299:   DM_DA          *da = (DM_DA*)dm->data;

1302:   if (da->bx != DM_BOUNDARY_PERIODIC) {
1303:     DMDAGetInfo(dm,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1304:     q    = (q-1)/(n-1);  /* number of spectral elements */
1305:     h    = 2.0/q;
1306:     DMDAGetCorners(dm,&xs,NULL,NULL,&xn,NULL,NULL);
1307:     xs   = xs/(n-1);
1308:     xn   = xn/(n-1);
1309:     DMDASetUniformCoordinates(dm,-1.,1.,0.,0.,0.,0.);
1310:     DMGetCoordinates(dm,&x);
1311:     DMDAVecGetArray(dm,x,&xx);

1313:     /* loop over local spectral elements */
1314:     for (j=xs; j<xs+xn; j++) {
1315:       /*
1316:        Except for the first process, each process starts on the second GLL point of the first element on that process
1317:        */
1318:       for (i= (j == xs && xs > 0)? 1 : 0; i<n; i++) {
1319:         xx[j*(n-1) + i] = -1.0 + h*j + h*(gll->nodes[i]+1.0)/2.;
1320:       }
1321:     }
1322:     DMDAVecRestoreArray(dm,x,&xx);
1323:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for periodic");
1324:   return(0);
1325: }

1327: /*@

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

1331:    Collective on DM

1333:    Input Parameters:
1334: +   da - the DMDA object
1335: -   gll - the GLL object

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

1341:    Level: advanced

1343: .seealso:   DMDACreate(), PetscGLLCreate(), DMGetCoordinates()
1344: @*/
1345: PetscErrorCode DMDASetGLLCoordinates(DM da,PetscGLL *gll)
1346: {

1350:   if (da->dim == 1) {
1351:     DMDASetGLLCoordinates_1d(da,gll);
1352:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for 2 or 3d");
1353:   return(0);
1354: }