Actual source code: da.c

  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;


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

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

 43:   Logically Collective on da

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

 51:   Level: intermediate

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

 64:   dd->m = m;
 65:   dd->n = n;
 66:   dd->p = p;
 67:   if (da->dim == 2) {
 68:     PetscMPIInt size;
 69:     MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
 70:     if ((dd->m > 0) && (dd->n < 0)) {
 71:       dd->n = size/dd->m;
 73:     }
 74:     if ((dd->n > 0) && (dd->m < 0)) {
 75:       dd->m = size/dd->n;
 77:     }
 78:   }
 79:   return 0;
 80: }

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

 85:   Not collective

 87:   Input Parameters:
 88: + da    - The DMDA
 89: - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC

 91:   Level: intermediate

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

104:   dd->bx = bx;
105:   dd->by = by;
106:   dd->bz = bz;
107:   return 0;
108: }

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

113:   Not collective

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

119:   Level: intermediate

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

130:   dd->w  = dof;
131:   da->bs = dof;
132:   return 0;
133: }

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

138:   Not collective

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

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

146:   Level: intermediate

148: .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
149: @*/
150: PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
151: {
152:   DM_DA *dd = (DM_DA *) da->data;

156:   *dof = dd->w;
157:   return 0;
158: }

160: /*@
161:   DMDAGetOverlap - Gets the size of the per-processor overlap.

163:   Not collective

165:   Input Parameter:
166: . da  - The DMDA

168:   Output Parameters:
169: + x   - Overlap in the x direction
170: . y   - Overlap in the y direction
171: - z   - Overlap in the z direction

173:   Level: intermediate

175: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
176: @*/
177: PetscErrorCode  DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
178: {
179:   DM_DA *dd = (DM_DA*)da->data;

182:   if (x) *x = dd->xol;
183:   if (y) *y = dd->yol;
184:   if (z) *z = dd->zol;
185:   return 0;
186: }

188: /*@
189:   DMDASetOverlap - Sets the size of the per-processor overlap.

191:   Not collective

193:   Input Parameters:
194: + da  - The DMDA
195: . x   - Overlap in the x direction
196: . y   - Overlap in the y direction
197: - z   - Overlap in the z direction

199:   Level: intermediate

201: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
202: @*/
203: PetscErrorCode  DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
204: {
205:   DM_DA *dd = (DM_DA*)da->data;

211:   dd->xol = x;
212:   dd->yol = y;
213:   dd->zol = z;
214:   return 0;
215: }

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

220:   Not collective

222:   Input Parameters:
223: . da  - The DMDA

225:   Output Parameters:
226: . Nsub   - Number of local subdomains created upon decomposition

228:   Level: intermediate

230: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
231: @*/
232: PetscErrorCode  DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
233: {
234:   DM_DA *dd = (DM_DA*)da->data;

237:   if (Nsub) *Nsub = dd->Nsub;
238:   return 0;
239: }

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

244:   Not collective

246:   Input Parameters:
247: + da  - The DMDA
248: - Nsub - The number of local subdomains requested

250:   Level: intermediate

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

260:   dd->Nsub = Nsub;
261:   return 0;
262: }

264: /*@
265:   DMDASetOffset - Sets the index offset of the DA.

267:   Collective on DA

269:   Input Parameters:
270: + da  - The DMDA
271: . xo  - The offset in the x direction
272: . yo  - The offset in the y direction
273: - zo  - The offset in the z direction

275:   Level: intermediate

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

281: .seealso: DMDAGetOffset(), DMDAVecGetArray()
282: @*/
283: PetscErrorCode  DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
284: {
285:   DM_DA          *dd = (DM_DA*)da->data;

294:   dd->xo = xo;
295:   dd->yo = yo;
296:   dd->zo = zo;
297:   dd->Mo = Mo;
298:   dd->No = No;
299:   dd->Po = Po;

301:   if (da->coordinateDM) {
302:     DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
303:   }
304:   return 0;
305: }

307: /*@
308:   DMDAGetOffset - Gets the index offset of the DA.

310:   Not collective

312:   Input Parameter:
313: . da  - The DMDA

315:   Output Parameters:
316: + xo  - The offset in the x direction
317: . yo  - The offset in the y direction
318: . zo  - The offset in the z direction
319: . Mo  - The global size in the x direction
320: . No  - The global size in the y direction
321: - Po  - The global size in the z direction

323:   Level: intermediate

325: .seealso: DMDASetOffset(), DMDAVecGetArray()
326: @*/
327: PetscErrorCode  DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
328: {
329:   DM_DA *dd = (DM_DA*)da->data;

332:   if (xo) *xo = dd->xo;
333:   if (yo) *yo = dd->yo;
334:   if (zo) *zo = dd->zo;
335:   if (Mo) *Mo = dd->Mo;
336:   if (No) *No = dd->No;
337:   if (Po) *Po = dd->Po;
338:   return 0;
339: }

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

344:   Not collective

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

349:   Output Parameters:
350: + xs  - The start of the region in x
351: . ys  - The start of the region in y
352: . zs  - The start of the region in z
353: . xs  - The size of the region in x
354: . ys  - The size of the region in y
355: - zs  - The size of the region in z

357:   Level: intermediate

359: .seealso: DMDAGetOffset(), DMDAVecGetArray()
360: @*/
361: PetscErrorCode  DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
362: {
363:   DM_DA          *dd = (DM_DA*)da->data;

366:   if (xs) *xs = dd->nonxs;
367:   if (ys) *ys = dd->nonys;
368:   if (zs) *zs = dd->nonzs;
369:   if (xm) *xm = dd->nonxm;
370:   if (ym) *ym = dd->nonym;
371:   if (zm) *zm = dd->nonzm;
372:   return 0;
373: }

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

378:   Collective on DA

380:   Input Parameters:
381: + da  - The DMDA
382: . xs  - The start of the region in x
383: . ys  - The start of the region in y
384: . zs  - The start of the region in z
385: . xs  - The size of the region in x
386: . ys  - The size of the region in y
387: - zs  - The size of the region in z

389:   Level: intermediate

391: .seealso: DMDAGetOffset(), DMDAVecGetArray()
392: @*/
393: PetscErrorCode  DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
394: {
395:   DM_DA          *dd = (DM_DA*)da->data;

404:   dd->nonxs = xs;
405:   dd->nonys = ys;
406:   dd->nonzs = zs;
407:   dd->nonxm = xm;
408:   dd->nonym = ym;
409:   dd->nonzm = zm;

411:   return 0;
412: }

414: /*@
415:   DMDASetStencilType - Sets the type of the communication stencil

417:   Logically Collective on da

419:   Input Parameters:
420: + da    - The DMDA
421: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.

423:   Level: intermediate

425: .seealso: DMDACreate(), DMDestroy(), DMDA
426: @*/
427: PetscErrorCode  DMDASetStencilType(DM da, DMDAStencilType stype)
428: {
429:   DM_DA *dd = (DM_DA*)da->data;

434:   dd->stencil_type = stype;
435:   return 0;
436: }

438: /*@
439:   DMDAGetStencilType - Gets the type of the communication stencil

441:   Not collective

443:   Input Parameter:
444: . da    - The DMDA

446:   Output Parameter:
447: . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.

449:   Level: intermediate

451: .seealso: DMDACreate(), DMDestroy(), DMDA
452: @*/
453: PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
454: {
455:   DM_DA *dd = (DM_DA*)da->data;

459:   *stype = dd->stencil_type;
460:   return 0;
461: }

463: /*@
464:   DMDASetStencilWidth - Sets the width of the communication stencil

466:   Logically Collective on da

468:   Input Parameters:
469: + da    - The DMDA
470: - width - The stencil width

472:   Level: intermediate

474: .seealso: DMDACreate(), DMDestroy(), DMDA
475: @*/
476: PetscErrorCode  DMDASetStencilWidth(DM da, PetscInt width)
477: {
478:   DM_DA *dd = (DM_DA*)da->data;

483:   dd->s = width;
484:   return 0;
485: }

487: /*@
488:   DMDAGetStencilWidth - Gets the width of the communication stencil

490:   Not collective

492:   Input Parameter:
493: . da    - The DMDA

495:   Output Parameter:
496: . width - The stencil width

498:   Level: intermediate

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

508:   *width = dd->s;
509:   return 0;
510: }

512: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
513: {
514:   PetscInt i,sum;

517:   for (i=sum=0; i<m; i++) sum += lx[i];
519:   return 0;
520: }

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

525:   Logically Collective on da

527:   Input Parameters:
528: + da - The DMDA
529: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
530: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
531: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.

533:   Level: intermediate

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

537: .seealso: DMDACreate(), DMDestroy(), DMDA
538: @*/
539: PetscErrorCode  DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
540: {
541:   DM_DA          *dd = (DM_DA*)da->data;

545:   if (lx) {
547:     DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
548:     if (!dd->lx) {
549:       PetscMalloc1(dd->m, &dd->lx);
550:     }
551:     PetscArraycpy(dd->lx, lx, dd->m);
552:   }
553:   if (ly) {
555:     DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
556:     if (!dd->ly) {
557:       PetscMalloc1(dd->n, &dd->ly);
558:     }
559:     PetscArraycpy(dd->ly, ly, dd->n);
560:   }
561:   if (lz) {
563:     DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
564:     if (!dd->lz) {
565:       PetscMalloc1(dd->p, &dd->lz);
566:     }
567:     PetscArraycpy(dd->lz, lz, dd->p);
568:   }
569:   return 0;
570: }

572: /*@
573:        DMDASetInterpolationType - Sets the type of interpolation that will be
574:           returned by DMCreateInterpolation()

576:    Logically Collective on da

578:    Input Parameters:
579: +  da - initial distributed array
580: -  ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms

582:    Level: intermediate

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

587: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
588: @*/
589: PetscErrorCode  DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
590: {
591:   DM_DA *dd = (DM_DA*)da->data;

595:   dd->interptype = ctype;
596:   return 0;
597: }

599: /*@
600:        DMDAGetInterpolationType - Gets the type of interpolation that will be
601:           used by DMCreateInterpolation()

603:    Not Collective

605:    Input Parameter:
606: .  da - distributed array

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

611:    Level: intermediate

613: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
614: @*/
615: PetscErrorCode  DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
616: {
617:   DM_DA *dd = (DM_DA*)da->data;

621:   *ctype = dd->interptype;
622:   return 0;
623: }

625: /*@C
626:       DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
627:         processes neighbors.

629:     Not Collective

631:    Input Parameter:
632: .     da - the DMDA object

634:    Output Parameters:
635: .     ranks - the neighbors ranks, stored with the x index increasing most rapidly.
636:               this process itself is in the list

638:    Notes:
639:     In 2d the array is of length 9, in 3d of length 27
640:           Not supported in 1d
641:           Do not free the array, it is freed when the DMDA is destroyed.

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

646:    Level: intermediate

648: @*/
649: PetscErrorCode  DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
650: {
651:   DM_DA *dd = (DM_DA*)da->data;

654:   *ranks = dd->neighbors;
655:   return 0;
656: }

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

661:     Not Collective

663:    Input Parameter:
664: .     da - the DMDA object

666:    Output Parameters:
667: +     lx - ownership along x direction (optional)
668: .     ly - ownership along y direction (optional)
669: -     lz - ownership along z direction (optional)

671:    Level: intermediate

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

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

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

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

683: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
684: @*/
685: PetscErrorCode  DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
686: {
687:   DM_DA *dd = (DM_DA*)da->data;

690:   if (lx) *lx = dd->lx;
691:   if (ly) *ly = dd->ly;
692:   if (lz) *lz = dd->lz;
693:   return 0;
694: }

696: /*@
697:      DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined

699:     Logically Collective on da

701:   Input Parameters:
702: +    da - the DMDA object
703: .    refine_x - ratio of fine grid to coarse in x direction (2 by default)
704: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
705: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

707:   Options Database:
708: +  -da_refine_x refine_x - refinement ratio in x direction
709: .  -da_refine_y rafine_y - refinement ratio in y direction
710: .  -da_refine_z refine_z - refinement ratio in z direction
711: -  -da_refine <n> - refine the DMDA object n times when it is created.

713:   Level: intermediate

715:     Notes:
716:     Pass PETSC_IGNORE to leave a value unchanged

718: .seealso: DMRefine(), DMDAGetRefinementFactor()
719: @*/
720: PetscErrorCode  DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
721: {
722:   DM_DA *dd = (DM_DA*)da->data;


729:   if (refine_x > 0) dd->refine_x = refine_x;
730:   if (refine_y > 0) dd->refine_y = refine_y;
731:   if (refine_z > 0) dd->refine_z = refine_z;
732:   return 0;
733: }

735: /*@C
736:      DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined

738:     Not Collective

740:   Input Parameter:
741: .    da - the DMDA object

743:   Output Parameters:
744: +    refine_x - ratio of fine grid to coarse in x direction (2 by default)
745: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
746: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

748:   Level: intermediate

750:     Notes:
751:     Pass NULL for values you do not need

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

760:   if (refine_x) *refine_x = dd->refine_x;
761:   if (refine_y) *refine_y = dd->refine_y;
762:   if (refine_z) *refine_z = dd->refine_z;
763:   return 0;
764: }

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

769:     Logically Collective on da

771:   Input Parameters:
772: +    da - the DMDA object
773: -    f - the function that allocates the matrix for that specific DMDA

775:   Level: developer

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

781:    Not supported from Fortran

783: .seealso: DMCreateMatrix(), DMDASetBlockFills()
784: @*/
785: PetscErrorCode  DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
786: {
788:   da->ops->creatematrix = f;
789:   return 0;
790: }

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

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

803:   if (ratio == 1) {
804:     PetscArraycpy(lf,lc,m);
805:     return 0;
806:   }
807:   for (i=0; i<m; i++) totalc += lc[i];
808:   remaining = (!periodic) + ratio * (totalc - (!periodic));
809:   for (i=0; i<m; i++) {
810:     PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
811:     if (i == m-1) lf[i] = want;
812:     else {
813:       const PetscInt nextc = startc + lc[i];
814:       /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
815:        * coarse stencil width of the first coarse node in the next subdomain. */
816:       while ((startf+want)/ratio < nextc - stencil_width) want++;
817:       /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
818:        * coarse stencil width of the last coarse node in the current subdomain. */
819:       while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
820:       /* Make sure all constraints are satisfied */
821:       if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
822:           || ((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");
823:     }
824:     lf[i]      = want;
825:     startc    += lc[i];
826:     startf    += lf[i];
827:     remaining -= lf[i];
828:   }
829:   return 0;
830: }

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

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

843:   if (ratio == 1) {
844:     PetscArraycpy(lc,lf,m);
845:     return 0;
846:   }
847:   for (i=0,totalf=0; i<m; i++) totalf += lf[i];
848:   remaining = (!periodic) + (totalf - (!periodic)) / ratio;
849:   for (i=0,startc=0,startf=0; i<m; i++) {
850:     PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
851:     if (i == m-1) lc[i] = want;
852:     else {
853:       const PetscInt nextf = startf+lf[i];
854:       /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
855:        * node is within one stencil width. */
856:       while (nextf/ratio < startc+want-stencil_width) want--;
857:       /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
858:        * fine node is within one stencil width. */
859:       while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
860:       if (want < 0 || want > remaining
861:           || (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");
862:     }
863:     lc[i]      = want;
864:     startc    += lc[i];
865:     startf    += lf[i];
866:     remaining -= lc[i];
867:   }
868:   return 0;
869: }

871: PetscErrorCode  DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
872: {
873:   PetscInt       M,N,P,i,dim;
874:   DM             da2;
875:   DM_DA          *dd = (DM_DA*)da->data,*dd2;


880:   DMGetDimension(da, &dim);
881:   if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
882:     M = dd->refine_x*dd->M;
883:   } else {
884:     M = 1 + dd->refine_x*(dd->M - 1);
885:   }
886:   if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
887:     if (dim > 1) {
888:       N = dd->refine_y*dd->N;
889:     } else {
890:       N = 1;
891:     }
892:   } else {
893:     N = 1 + dd->refine_y*(dd->N - 1);
894:   }
895:   if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
896:     if (dim > 2) {
897:       P = dd->refine_z*dd->P;
898:     } else {
899:       P = 1;
900:     }
901:   } else {
902:     P = 1 + dd->refine_z*(dd->P - 1);
903:   }
904:   DMDACreate(PetscObjectComm((PetscObject)da),&da2);
905:   DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
906:   DMSetDimension(da2,dim);
907:   DMDASetSizes(da2,M,N,P);
908:   DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
909:   DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
910:   DMDASetDof(da2,dd->w);
911:   DMDASetStencilType(da2,dd->stencil_type);
912:   DMDASetStencilWidth(da2,dd->s);
913:   if (dim == 3) {
914:     PetscInt *lx,*ly,*lz;
915:     PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
916:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
917:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
918:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
919:     DMDASetOwnershipRanges(da2,lx,ly,lz);
920:     PetscFree3(lx,ly,lz);
921:   } else if (dim == 2) {
922:     PetscInt *lx,*ly;
923:     PetscMalloc2(dd->m,&lx,dd->n,&ly);
924:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
925:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
926:     DMDASetOwnershipRanges(da2,lx,ly,NULL);
927:     PetscFree2(lx,ly);
928:   } else if (dim == 1) {
929:     PetscInt *lx;
930:     PetscMalloc1(dd->m,&lx);
931:     DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
932:     DMDASetOwnershipRanges(da2,lx,NULL,NULL);
933:     PetscFree(lx);
934:   }
935:   dd2 = (DM_DA*)da2->data;

937:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
938:   da2->ops->creatematrix = da->ops->creatematrix;
939:   /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
940:   da2->ops->getcoloring = da->ops->getcoloring;
941:   dd2->interptype       = dd->interptype;

943:   /* copy fill information if given */
944:   if (dd->dfill) {
945:     PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
946:     PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);
947:   }
948:   if (dd->ofill) {
949:     PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
950:     PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);
951:   }
952:   /* copy the refine information */
953:   dd2->coarsen_x = dd2->refine_x = dd->refine_x;
954:   dd2->coarsen_y = dd2->refine_y = dd->refine_y;
955:   dd2->coarsen_z = dd2->refine_z = dd->refine_z;

957:   if (dd->refine_z_hier) {
958:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) {
959:       dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
960:     }
961:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) {
962:       dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
963:     }
964:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
965:     PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
966:     PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);
967:   }
968:   if (dd->refine_y_hier) {
969:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) {
970:       dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
971:     }
972:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) {
973:       dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
974:     }
975:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
976:     PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
977:     PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);
978:   }
979:   if (dd->refine_x_hier) {
980:     if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) {
981:       dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
982:     }
983:     if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) {
984:       dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
985:     }
986:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
987:     PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
988:     PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);
989:   }

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

994:   dd2->lf = dd->lf;
995:   dd2->lj = dd->lj;

997:   da2->leveldown = da->leveldown;
998:   da2->levelup   = da->levelup + 1;

1000:   DMSetUp(da2);

1002:   /* interpolate coordinates if they are set on the coarse grid */
1003:   if (da->coordinates) {
1004:     DM  cdaf,cdac;
1005:     Vec coordsc,coordsf;
1006:     Mat II;

1008:     DMGetCoordinateDM(da,&cdac);
1009:     DMGetCoordinates(da,&coordsc);
1010:     DMGetCoordinateDM(da2,&cdaf);
1011:     /* force creation of the coordinate vector */
1012:     DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1013:     DMGetCoordinates(da2,&coordsf);
1014:     DMCreateInterpolation(cdac,cdaf,&II,NULL);
1015:     MatInterpolate(II,coordsc,coordsf);
1016:     MatDestroy(&II);
1017:   }

1019:   for (i=0; i<da->bs; i++) {
1020:     const char *fieldname;
1021:     DMDAGetFieldName(da,i,&fieldname);
1022:     DMDASetFieldName(da2,i,fieldname);
1023:   }

1025:   *daref = da2;
1026:   return 0;
1027: }

1029: PetscErrorCode  DMCoarsen_DA(DM dmf, MPI_Comm comm,DM *dmc)
1030: {
1031:   PetscInt       M,N,P,i,dim;
1032:   DM             dmc2;
1033:   DM_DA          *dd = (DM_DA*)dmf->data,*dd2;


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

1095:   /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1096:   /* dmc2->ops->createinterpolation = dmf->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1097:   dmc2->ops->creatematrix = dmf->ops->creatematrix;
1098:   dmc2->ops->getcoloring  = dmf->ops->getcoloring;
1099:   dd2->interptype         = dd->interptype;

1101:   /* copy fill information if given */
1102:   if (dd->dfill) {
1103:     PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1104:     PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);
1105:   }
1106:   if (dd->ofill) {
1107:     PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1108:     PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);
1109:   }
1110:   /* copy the refine information */
1111:   dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1112:   dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1113:   dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;

1115:   if (dd->refine_z_hier) {
1116:     if (dmf->levelup - dmf->leveldown -1 > -1 && dmf->levelup - dmf->leveldown - 1< dd->refine_z_hier_n) {
1117:       dd2->refine_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 1];
1118:     }
1119:     if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_z_hier_n) {
1120:       dd2->coarsen_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 2];
1121:     }
1122:     dd2->refine_z_hier_n = dd->refine_z_hier_n;
1123:     PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);
1124:     PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);
1125:   }
1126:   if (dd->refine_y_hier) {
1127:     if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1< dd->refine_y_hier_n) {
1128:       dd2->refine_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 1];
1129:     }
1130:     if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_y_hier_n) {
1131:       dd2->coarsen_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 2];
1132:     }
1133:     dd2->refine_y_hier_n = dd->refine_y_hier_n;
1134:     PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);
1135:     PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);
1136:   }
1137:   if (dd->refine_x_hier) {
1138:     if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_x_hier_n) {
1139:       dd2->refine_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 1];
1140:     }
1141:     if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_x_hier_n) {
1142:       dd2->coarsen_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 2];
1143:     }
1144:     dd2->refine_x_hier_n = dd->refine_x_hier_n;
1145:     PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);
1146:     PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);
1147:   }

1149:   /* copy vector type information */
1150:   DMSetVecType(dmc2,dmf->vectype);

1152:   dd2->lf = dd->lf;
1153:   dd2->lj = dd->lj;

1155:   dmc2->leveldown = dmf->leveldown + 1;
1156:   dmc2->levelup   = dmf->levelup;

1158:   DMSetUp(dmc2);

1160:   /* inject coordinates if they are set on the fine grid */
1161:   if (dmf->coordinates) {
1162:     DM         cdaf,cdac;
1163:     Vec        coordsc,coordsf;
1164:     Mat        inject;
1165:     VecScatter vscat;

1167:     DMGetCoordinateDM(dmf,&cdaf);
1168:     DMGetCoordinates(dmf,&coordsf);
1169:     DMGetCoordinateDM(dmc2,&cdac);
1170:     /* force creation of the coordinate vector */
1171:     DMDASetUniformCoordinates(dmc2,0.0,1.0,0.0,1.0,0.0,1.0);
1172:     DMGetCoordinates(dmc2,&coordsc);

1174:     DMCreateInjection(cdac,cdaf,&inject);
1175:     MatScatterGetVecScatter(inject,&vscat);
1176:     VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1177:     VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1178:     MatDestroy(&inject);
1179:   }

1181:   for (i=0; i<dmf->bs; i++) {
1182:     const char *fieldname;
1183:     DMDAGetFieldName(dmf,i,&fieldname);
1184:     DMDASetFieldName(dmc2,i,fieldname);
1185:   }

1187:   *dmc = dmc2;
1188:   return 0;
1189: }

1191: PetscErrorCode  DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1192: {
1193:   PetscInt       i,n,*refx,*refy,*refz;

1197:   if (nlevels == 0) return 0;

1200:   /* Get refinement factors, defaults taken from the coarse DMDA */
1201:   PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1202:   for (i=0; i<nlevels; i++) {
1203:     DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1204:   }
1205:   n    = nlevels;
1206:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1207:   n    = nlevels;
1208:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1209:   n    = nlevels;
1210:   PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);

1212:   DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1213:   DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1214:   for (i=1; i<nlevels; i++) {
1215:     DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1216:     DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1217:   }
1218:   PetscFree3(refx,refy,refz);
1219:   return 0;
1220: }

1222: PetscErrorCode  DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1223: {
1224:   PetscInt       i;

1228:   if (nlevels == 0) return 0;
1230:   DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1231:   for (i=1; i<nlevels; i++) {
1232:     DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1233:   }
1234:   return 0;
1235: }

1237: PetscErrorCode DMDASetGLLCoordinates_1d(DM dm,PetscInt n,PetscReal *nodes)
1238: {
1239:   PetscInt       i,j,xs,xn,q;
1240:   PetscScalar    *xx;
1241:   PetscReal      h;
1242:   Vec            x;
1243:   DM_DA          *da = (DM_DA*)dm->data;

1245:   if (da->bx != DM_BOUNDARY_PERIODIC) {
1246:     DMDAGetInfo(dm,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1247:     q    = (q-1)/(n-1);  /* number of spectral elements */
1248:     h    = 2.0/q;
1249:     DMDAGetCorners(dm,&xs,NULL,NULL,&xn,NULL,NULL);
1250:     xs   = xs/(n-1);
1251:     xn   = xn/(n-1);
1252:     DMDASetUniformCoordinates(dm,-1.,1.,0.,0.,0.,0.);
1253:     DMGetCoordinates(dm,&x);
1254:     DMDAVecGetArray(dm,x,&xx);

1256:     /* loop over local spectral elements */
1257:     for (j=xs; j<xs+xn; j++) {
1258:       /*
1259:        Except for the first process, each process starts on the second GLL point of the first element on that process
1260:        */
1261:       for (i= (j == xs && xs > 0)? 1 : 0; i<n; i++) {
1262:         xx[j*(n-1) + i] = -1.0 + h*j + h*(nodes[i]+1.0)/2.;
1263:       }
1264:     }
1265:     DMDAVecRestoreArray(dm,x,&xx);
1266:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for periodic");
1267:   return 0;
1268: }

1270: /*@

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

1274:    Collective on da

1276:    Input Parameters:
1277: +   da - the DMDA object
1278: -   n - the number of GLL nodes
1279: -   nodes - the GLL nodes

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

1286:    Level: advanced

1288: .seealso:   DMDACreate(), PetscDTGaussLobattoLegendreQuadrature(), DMGetCoordinates()
1289: @*/
1290: PetscErrorCode DMDASetGLLCoordinates(DM da,PetscInt n,PetscReal *nodes)
1291: {
1292:   if (da->dim == 1) {
1293:     DMDASetGLLCoordinates_1d(da,n,nodes);
1294:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for 2 or 3d");
1295:   return 0;
1296: }

1298: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM da1,DM dm2,PetscBool *compatible,PetscBool *set)
1299: {
1300:   DM_DA          *dd1 = (DM_DA*)da1->data,*dd2;
1301:   DM             da2;
1302:   DMType         dmtype2;
1303:   PetscBool      isda,compatibleLocal;
1304:   PetscInt       i;

1307:   DMGetType(dm2,&dmtype2);
1308:   PetscStrcmp(dmtype2,DMDA,&isda);
1309:   if (isda) {
1310:     da2 = dm2;
1311:     dd2 = (DM_DA*)da2->data;
1313:     compatibleLocal = (PetscBool)(da1->dim == da2->dim);
1314:     if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
1315:     /*                                                                           Global size              ranks               Boundary type */
1316:     if (compatibleLocal)                 compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
1317:     if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
1318:     if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
1319:     if (compatibleLocal) {
1320:       for (i=0; i<dd1->m; ++i) {
1321:         compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i]));           /* Local size     */
1322:       }
1323:     }
1324:     if (compatibleLocal && da1->dim > 1) {
1325:       for (i=0; i<dd1->n; ++i) {
1326:         compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
1327:       }
1328:     }
1329:     if (compatibleLocal && da1->dim > 2) {
1330:       for (i=0; i<dd1->p; ++i) {
1331:         compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
1332:       }
1333:     }
1334:     *compatible = compatibleLocal;
1335:     *set = PETSC_TRUE;
1336:   } else {
1337:     /* Decline to determine compatibility with other DM types */
1338:     *set = PETSC_FALSE;
1339:   }
1340:   return 0;
1341: }