Actual source code: dacreate.c


  2: #include <petsc/private/dmdaimpl.h>

  4: PetscErrorCode DMSetFromOptions_DA(DM da, PetscOptionItems *PetscOptionsObject)
  5: {
  6:   DM_DA    *dd     = (DM_DA *)da->data;
  7:   PetscInt  refine = 0, dim = da->dim, maxnlevels = 100, refx[100], refy[100], refz[100], n, i;
  8:   PetscBool flg;


 14:   PetscOptionsHeadBegin(PetscOptionsObject, "DMDA Options");
 15:   PetscOptionsBoundedInt("-da_grid_x", "Number of grid points in x direction", "DMDASetSizes", dd->M, &dd->M, NULL, 1);
 16:   if (dim > 1) PetscOptionsBoundedInt("-da_grid_y", "Number of grid points in y direction", "DMDASetSizes", dd->N, &dd->N, NULL, 1);
 17:   if (dim > 2) PetscOptionsBoundedInt("-da_grid_z", "Number of grid points in z direction", "DMDASetSizes", dd->P, &dd->P, NULL, 1);

 19:   PetscOptionsBoundedInt("-da_overlap", "Decomposition overlap in all directions", "DMDASetOverlap", dd->xol, &dd->xol, &flg, 0);
 20:   if (flg) DMDASetOverlap(da, dd->xol, dd->xol, dd->xol);
 21:   PetscOptionsBoundedInt("-da_overlap_x", "Decomposition overlap in x direction", "DMDASetOverlap", dd->xol, &dd->xol, NULL, 0);
 22:   if (dim > 1) PetscOptionsBoundedInt("-da_overlap_y", "Decomposition overlap in y direction", "DMDASetOverlap", dd->yol, &dd->yol, NULL, 0);
 23:   if (dim > 2) PetscOptionsBoundedInt("-da_overlap_z", "Decomposition overlap in z direction", "DMDASetOverlap", dd->zol, &dd->zol, NULL, 0);

 25:   PetscOptionsBoundedInt("-da_local_subdomains", "", "DMDASetNumLocalSubdomains", dd->Nsub, &dd->Nsub, &flg, PETSC_DECIDE);
 26:   if (flg) DMDASetNumLocalSubDomains(da, dd->Nsub);

 28:   /* Handle DMDA parallel distribution */
 29:   PetscOptionsBoundedInt("-da_processors_x", "Number of processors in x direction", "DMDASetNumProcs", dd->m, &dd->m, NULL, PETSC_DECIDE);
 30:   if (dim > 1) PetscOptionsBoundedInt("-da_processors_y", "Number of processors in y direction", "DMDASetNumProcs", dd->n, &dd->n, NULL, PETSC_DECIDE);
 31:   if (dim > 2) PetscOptionsBoundedInt("-da_processors_z", "Number of processors in z direction", "DMDASetNumProcs", dd->p, &dd->p, NULL, PETSC_DECIDE);
 32:   /* Handle DMDA refinement */
 33:   PetscOptionsBoundedInt("-da_refine_x", "Refinement ratio in x direction", "DMDASetRefinementFactor", dd->refine_x, &dd->refine_x, NULL, 1);
 34:   if (dim > 1) PetscOptionsBoundedInt("-da_refine_y", "Refinement ratio in y direction", "DMDASetRefinementFactor", dd->refine_y, &dd->refine_y, NULL, 1);
 35:   if (dim > 2) PetscOptionsBoundedInt("-da_refine_z", "Refinement ratio in z direction", "DMDASetRefinementFactor", dd->refine_z, &dd->refine_z, NULL, 1);
 36:   dd->coarsen_x = dd->refine_x;
 37:   dd->coarsen_y = dd->refine_y;
 38:   dd->coarsen_z = dd->refine_z;

 40:   /* Get refinement factors, defaults taken from the coarse DMDA */
 41:   DMDAGetRefinementFactor(da, &refx[0], &refy[0], &refz[0]);
 42:   for (i = 1; i < maxnlevels; i++) {
 43:     refx[i] = refx[0];
 44:     refy[i] = refy[0];
 45:     refz[i] = refz[0];
 46:   }
 47:   n = maxnlevels;
 48:   PetscOptionsIntArray("-da_refine_hierarchy_x", "Refinement factor for each level", "None", refx, &n, &flg);
 49:   if (flg) {
 50:     dd->refine_x        = refx[0];
 51:     dd->refine_x_hier_n = n;
 52:     PetscMalloc1(n, &dd->refine_x_hier);
 53:     PetscArraycpy(dd->refine_x_hier, refx, n);
 54:   }
 55:   if (dim > 1) {
 56:     n = maxnlevels;
 57:     PetscOptionsIntArray("-da_refine_hierarchy_y", "Refinement factor for each level", "None", refy, &n, &flg);
 58:     if (flg) {
 59:       dd->refine_y        = refy[0];
 60:       dd->refine_y_hier_n = n;
 61:       PetscMalloc1(n, &dd->refine_y_hier);
 62:       PetscArraycpy(dd->refine_y_hier, refy, n);
 63:     }
 64:   }
 65:   if (dim > 2) {
 66:     n = maxnlevels;
 67:     PetscOptionsIntArray("-da_refine_hierarchy_z", "Refinement factor for each level", "None", refz, &n, &flg);
 68:     if (flg) {
 69:       dd->refine_z        = refz[0];
 70:       dd->refine_z_hier_n = n;
 71:       PetscMalloc1(n, &dd->refine_z_hier);
 72:       PetscArraycpy(dd->refine_z_hier, refz, n);
 73:     }
 74:   }

 76:   PetscOptionsBoundedInt("-da_refine", "Uniformly refine DA one or more times", "None", refine, &refine, NULL, 0);
 77:   PetscOptionsHeadEnd();

 79:   while (refine--) {
 80:     if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
 81:       PetscIntMultError(dd->refine_x, dd->M, &dd->M);
 82:     } else {
 83:       PetscIntMultError(dd->refine_x, dd->M - 1, &dd->M);
 84:       dd->M += 1;
 85:     }
 86:     if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
 87:       PetscIntMultError(dd->refine_y, dd->N, &dd->N);
 88:     } else {
 89:       PetscIntMultError(dd->refine_y, dd->N - 1, &dd->N);
 90:       dd->N += 1;
 91:     }
 92:     if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
 93:       PetscIntMultError(dd->refine_z, dd->P, &dd->P);
 94:     } else {
 95:       PetscIntMultError(dd->refine_z, dd->P - 1, &dd->P);
 96:       dd->P += 1;
 97:     }
 98:     da->levelup++;
 99:     if (da->levelup - da->leveldown >= 0) {
100:       dd->refine_x = refx[da->levelup - da->leveldown];
101:       dd->refine_y = refy[da->levelup - da->leveldown];
102:       dd->refine_z = refz[da->levelup - da->leveldown];
103:     }
104:     if (da->levelup - da->leveldown >= 1) {
105:       dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
106:       dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
107:       dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
108:     }
109:   }
110:   return 0;
111: }

113: extern PetscErrorCode       DMCreateGlobalVector_DA(DM, Vec *);
114: extern PetscErrorCode       DMCreateLocalVector_DA(DM, Vec *);
115: extern PetscErrorCode       DMGlobalToLocalBegin_DA(DM, Vec, InsertMode, Vec);
116: extern PetscErrorCode       DMGlobalToLocalEnd_DA(DM, Vec, InsertMode, Vec);
117: extern PetscErrorCode       DMLocalToGlobalBegin_DA(DM, Vec, InsertMode, Vec);
118: extern PetscErrorCode       DMLocalToGlobalEnd_DA(DM, Vec, InsertMode, Vec);
119: extern PetscErrorCode       DMLocalToLocalBegin_DA(DM, Vec, InsertMode, Vec);
120: extern PetscErrorCode       DMLocalToLocalEnd_DA(DM, Vec, InsertMode, Vec);
121: extern PetscErrorCode       DMCreateInterpolation_DA(DM, DM, Mat *, Vec *);
122: extern PetscErrorCode       DMCreateColoring_DA(DM, ISColoringType, ISColoring *);
123: extern PetscErrorCode       DMCreateMatrix_DA(DM, Mat *);
124: extern PetscErrorCode       DMCreateCoordinateDM_DA(DM, DM *);
125: extern PetscErrorCode       DMRefine_DA(DM, MPI_Comm, DM *);
126: extern PetscErrorCode       DMCoarsen_DA(DM, MPI_Comm, DM *);
127: extern PetscErrorCode       DMRefineHierarchy_DA(DM, PetscInt, DM[]);
128: extern PetscErrorCode       DMCoarsenHierarchy_DA(DM, PetscInt, DM[]);
129: extern PetscErrorCode       DMCreateInjection_DA(DM, DM, Mat *);
130: extern PetscErrorCode       DMView_DA(DM, PetscViewer);
131: extern PetscErrorCode       DMSetUp_DA(DM);
132: extern PetscErrorCode       DMDestroy_DA(DM);
133: extern PetscErrorCode       DMCreateDomainDecomposition_DA(DM, PetscInt *, char ***, IS **, IS **, DM **);
134: extern PetscErrorCode       DMCreateDomainDecompositionScatters_DA(DM, PetscInt, DM *, VecScatter **, VecScatter **, VecScatter **);
135: PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM, DM, PetscBool *, PetscBool *);

137: PetscErrorCode DMLoad_DA(DM da, PetscViewer viewer)
138: {
139:   PetscInt        dim, m, n, p, dof, swidth;
140:   DMDAStencilType stencil;
141:   DMBoundaryType  bx, by, bz;
142:   PetscBool       coors;
143:   DM              dac;
144:   Vec             c;

146:   PetscViewerBinaryRead(viewer, &dim, 1, NULL, PETSC_INT);
147:   PetscViewerBinaryRead(viewer, &m, 1, NULL, PETSC_INT);
148:   PetscViewerBinaryRead(viewer, &n, 1, NULL, PETSC_INT);
149:   PetscViewerBinaryRead(viewer, &p, 1, NULL, PETSC_INT);
150:   PetscViewerBinaryRead(viewer, &dof, 1, NULL, PETSC_INT);
151:   PetscViewerBinaryRead(viewer, &swidth, 1, NULL, PETSC_INT);
152:   PetscViewerBinaryRead(viewer, &bx, 1, NULL, PETSC_ENUM);
153:   PetscViewerBinaryRead(viewer, &by, 1, NULL, PETSC_ENUM);
154:   PetscViewerBinaryRead(viewer, &bz, 1, NULL, PETSC_ENUM);
155:   PetscViewerBinaryRead(viewer, &stencil, 1, NULL, PETSC_ENUM);

157:   DMSetDimension(da, dim);
158:   DMDASetSizes(da, m, n, p);
159:   DMDASetBoundaryType(da, bx, by, bz);
160:   DMDASetDof(da, dof);
161:   DMDASetStencilType(da, stencil);
162:   DMDASetStencilWidth(da, swidth);
163:   DMSetUp(da);
164:   PetscViewerBinaryRead(viewer, &coors, 1, NULL, PETSC_ENUM);
165:   if (coors) {
166:     DMGetCoordinateDM(da, &dac);
167:     DMCreateGlobalVector(dac, &c);
168:     VecLoad(c, viewer);
169:     DMSetCoordinates(da, c);
170:     VecDestroy(&c);
171:   }
172:   return 0;
173: }

175: PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
176: {
177:   DM_DA *da = (DM_DA *)dm->data;

179:   if (subdm) {
180:     PetscSF sf;
181:     Vec     coords;
182:     void   *ctx;
183:     /* Cannot use DMClone since the dof stuff is mixed in. Ugh
184:     DMClone(dm, subdm); */
185:     DMCreate(PetscObjectComm((PetscObject)dm), subdm);
186:     DMGetPointSF(dm, &sf);
187:     DMSetPointSF(*subdm, sf);
188:     DMGetApplicationContext(dm, &ctx);
189:     DMSetApplicationContext(*subdm, ctx);
190:     DMGetCoordinatesLocal(dm, &coords);
191:     if (coords) {
192:       DMSetCoordinatesLocal(*subdm, coords);
193:     } else {
194:       DMGetCoordinates(dm, &coords);
195:       if (coords) DMSetCoordinates(*subdm, coords);
196:     }

198:     DMSetType(*subdm, DMDA);
199:     DMSetDimension(*subdm, dm->dim);
200:     DMDASetSizes(*subdm, da->M, da->N, da->P);
201:     DMDASetNumProcs(*subdm, da->m, da->n, da->p);
202:     DMDASetBoundaryType(*subdm, da->bx, da->by, da->bz);
203:     DMDASetDof(*subdm, numFields);
204:     DMDASetStencilType(*subdm, da->stencil_type);
205:     DMDASetStencilWidth(*subdm, da->s);
206:     DMDASetOwnershipRanges(*subdm, da->lx, da->ly, da->lz);
207:   }
208:   if (is) {
209:     PetscInt *indices, cnt = 0, dof = da->w, i, j;

211:     PetscMalloc1(da->Nlocal * numFields / dof, &indices);
212:     for (i = da->base / dof; i < (da->base + da->Nlocal) / dof; ++i) {
213:       for (j = 0; j < numFields; ++j) indices[cnt++] = dof * i + fields[j];
214:     }
216:     ISCreateGeneral(PetscObjectComm((PetscObject)dm), cnt, indices, PETSC_OWN_POINTER, is);
217:   }
218:   return 0;
219: }

221: PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
222: {
223:   PetscInt i;
224:   DM_DA   *dd  = (DM_DA *)dm->data;
225:   PetscInt dof = dd->w;

227:   if (len) *len = dof;
228:   if (islist) {
229:     Vec      v;
230:     PetscInt rstart, n;

232:     DMGetGlobalVector(dm, &v);
233:     VecGetOwnershipRange(v, &rstart, NULL);
234:     VecGetLocalSize(v, &n);
235:     DMRestoreGlobalVector(dm, &v);
236:     PetscMalloc1(dof, islist);
237:     for (i = 0; i < dof; i++) ISCreateStride(PetscObjectComm((PetscObject)dm), n / dof, rstart + i, dof, &(*islist)[i]);
238:   }
239:   if (namelist) {
240:     PetscMalloc1(dof, namelist);
241:     if (dd->fieldname) {
242:       for (i = 0; i < dof; i++) PetscStrallocpy(dd->fieldname[i], &(*namelist)[i]);
243:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Currently DMDA must have fieldnames");
244:   }
245:   if (dmlist) {
246:     DM da;

248:     DMDACreate(PetscObjectComm((PetscObject)dm), &da);
249:     DMSetDimension(da, dm->dim);
250:     DMDASetSizes(da, dd->M, dd->N, dd->P);
251:     DMDASetNumProcs(da, dd->m, dd->n, dd->p);
252:     DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);
253:     DMDASetDof(da, 1);
254:     DMDASetStencilType(da, dd->stencil_type);
255:     DMDASetStencilWidth(da, dd->s);
256:     DMSetUp(da);
257:     PetscMalloc1(dof, dmlist);
258:     for (i = 0; i < dof - 1; i++) PetscObjectReference((PetscObject)da);
259:     for (i = 0; i < dof; i++) (*dmlist)[i] = da;
260:   }
261:   return 0;
262: }

264: PetscErrorCode DMClone_DA(DM dm, DM *newdm)
265: {
266:   DM_DA *da = (DM_DA *)dm->data;

268:   DMSetType(*newdm, DMDA);
269:   DMSetDimension(*newdm, dm->dim);
270:   DMDASetSizes(*newdm, da->M, da->N, da->P);
271:   DMDASetNumProcs(*newdm, da->m, da->n, da->p);
272:   DMDASetBoundaryType(*newdm, da->bx, da->by, da->bz);
273:   DMDASetDof(*newdm, da->w);
274:   DMDASetStencilType(*newdm, da->stencil_type);
275:   DMDASetStencilWidth(*newdm, da->s);
276:   DMDASetOwnershipRanges(*newdm, da->lx, da->ly, da->lz);
277:   DMSetUp(*newdm);
278:   return 0;
279: }

281: static PetscErrorCode DMHasCreateInjection_DA(DM dm, PetscBool *flg)
282: {
283:   DM_DA *da = (DM_DA *)dm->data;

287:   *flg = da->interptype == DMDA_Q1 ? PETSC_TRUE : PETSC_FALSE;
288:   return 0;
289: }

291: static PetscErrorCode DMGetDimPoints_DA(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
292: {
293:   DMDAGetDepthStratum(dm, dim, pStart, pEnd);
294:   return 0;
295: }

297: static PetscErrorCode DMGetNeighbors_DA(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[])
298: {
299:   PetscInt        dim;
300:   DMDAStencilType st;

302:   DMDAGetNeighbors(dm, ranks);
303:   DMDAGetInfo(dm, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &st);

305:   switch (dim) {
306:   case 1:
307:     *nranks = 3;
308:     /* if (st == DMDA_STENCIL_STAR) *nranks = 3; */
309:     break;
310:   case 2:
311:     *nranks = 9;
312:     /* if (st == DMDA_STENCIL_STAR) *nranks = 5; */
313:     break;
314:   case 3:
315:     *nranks = 27;
316:     /* if (st == DMDA_STENCIL_STAR) *nranks = 7; */
317:     break;
318:   default:
319:     break;
320:   }
321:   return 0;
322: }

324: /*MC
325:    DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions.
326:          In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points.
327:          In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width.

329:          The vectors can be thought of as either cell centered or vertex centered on the mesh. But some variables cannot be cell centered and others
330:          vertex centered; see the documentation for DMSTAG, a similar DM implementation which supports these staggered grids.

332:   Level: intermediate

334: .seealso: `DMType`, `DMCOMPOSITE`, `DMSTAG`, `DMDACreate()`, `DMCreate()`, `DMSetType()`
335: M*/

337: extern PetscErrorCode       DMLocatePoints_DA_Regular(DM, Vec, DMPointLocationType, PetscSF);
338: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_DMDA(PetscObject, PetscViewer);

340: PETSC_EXTERN PetscErrorCode DMCreate_DA(DM da)
341: {
342:   DM_DA *dd;

345:   PetscNew(&dd);
346:   da->data = dd;

348:   da->dim        = -1;
349:   dd->interptype = DMDA_Q1;
350:   dd->refine_x   = 2;
351:   dd->refine_y   = 2;
352:   dd->refine_z   = 2;
353:   dd->coarsen_x  = 2;
354:   dd->coarsen_y  = 2;
355:   dd->coarsen_z  = 2;
356:   dd->fieldname  = NULL;
357:   dd->nlocal     = -1;
358:   dd->Nlocal     = -1;
359:   dd->M          = -1;
360:   dd->N          = -1;
361:   dd->P          = -1;
362:   dd->m          = -1;
363:   dd->n          = -1;
364:   dd->p          = -1;
365:   dd->w          = -1;
366:   dd->s          = -1;

368:   dd->xs = -1;
369:   dd->xe = -1;
370:   dd->ys = -1;
371:   dd->ye = -1;
372:   dd->zs = -1;
373:   dd->ze = -1;
374:   dd->Xs = -1;
375:   dd->Xe = -1;
376:   dd->Ys = -1;
377:   dd->Ye = -1;
378:   dd->Zs = -1;
379:   dd->Ze = -1;

381:   dd->Nsub = 1;
382:   dd->xol  = 0;
383:   dd->yol  = 0;
384:   dd->zol  = 0;
385:   dd->xo   = 0;
386:   dd->yo   = 0;
387:   dd->zo   = 0;
388:   dd->Mo   = -1;
389:   dd->No   = -1;
390:   dd->Po   = -1;

392:   dd->gtol = NULL;
393:   dd->ltol = NULL;
394:   dd->ao   = NULL;
395:   PetscStrallocpy(AOBASIC, (char **)&dd->aotype);
396:   dd->base         = -1;
397:   dd->bx           = DM_BOUNDARY_NONE;
398:   dd->by           = DM_BOUNDARY_NONE;
399:   dd->bz           = DM_BOUNDARY_NONE;
400:   dd->stencil_type = DMDA_STENCIL_BOX;
401:   dd->interptype   = DMDA_Q1;
402:   dd->lx           = NULL;
403:   dd->ly           = NULL;
404:   dd->lz           = NULL;

406:   dd->elementtype = DMDA_ELEMENT_Q1;

408:   da->ops->globaltolocalbegin        = DMGlobalToLocalBegin_DA;
409:   da->ops->globaltolocalend          = DMGlobalToLocalEnd_DA;
410:   da->ops->localtoglobalbegin        = DMLocalToGlobalBegin_DA;
411:   da->ops->localtoglobalend          = DMLocalToGlobalEnd_DA;
412:   da->ops->localtolocalbegin         = DMLocalToLocalBegin_DA;
413:   da->ops->localtolocalend           = DMLocalToLocalEnd_DA;
414:   da->ops->createglobalvector        = DMCreateGlobalVector_DA;
415:   da->ops->createlocalvector         = DMCreateLocalVector_DA;
416:   da->ops->createinterpolation       = DMCreateInterpolation_DA;
417:   da->ops->getcoloring               = DMCreateColoring_DA;
418:   da->ops->creatematrix              = DMCreateMatrix_DA;
419:   da->ops->refine                    = DMRefine_DA;
420:   da->ops->coarsen                   = DMCoarsen_DA;
421:   da->ops->refinehierarchy           = DMRefineHierarchy_DA;
422:   da->ops->coarsenhierarchy          = DMCoarsenHierarchy_DA;
423:   da->ops->createinjection           = DMCreateInjection_DA;
424:   da->ops->hascreateinjection        = DMHasCreateInjection_DA;
425:   da->ops->destroy                   = DMDestroy_DA;
426:   da->ops->view                      = NULL;
427:   da->ops->setfromoptions            = DMSetFromOptions_DA;
428:   da->ops->setup                     = DMSetUp_DA;
429:   da->ops->clone                     = DMClone_DA;
430:   da->ops->load                      = DMLoad_DA;
431:   da->ops->createcoordinatedm        = DMCreateCoordinateDM_DA;
432:   da->ops->createsubdm               = DMCreateSubDM_DA;
433:   da->ops->createfielddecomposition  = DMCreateFieldDecomposition_DA;
434:   da->ops->createdomaindecomposition = DMCreateDomainDecomposition_DA;
435:   da->ops->createddscatters          = DMCreateDomainDecompositionScatters_DA;
436:   da->ops->getdimpoints              = DMGetDimPoints_DA;
437:   da->ops->getneighbors              = DMGetNeighbors_DA;
438:   da->ops->locatepoints              = DMLocatePoints_DA_Regular;
439:   da->ops->getcompatibility          = DMGetCompatibility_DA;
440:   PetscObjectComposeFunction((PetscObject)da, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_DMDA);
441:   return 0;
442: }

444: /*@
445:   DMDACreate - Creates a DMDA object.

447:   Collective

449:   Input Parameter:
450: . comm - The communicator for the DMDA object

452:   Output Parameter:
453: . da  - The DMDA object

455:   Level: advanced

457:   Developers Note:
458:   Since there exists DMDACreate1/2/3d() should this routine even exist?

460: .seealso: `DMDASetSizes()`, `DMClone()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
461: @*/
462: PetscErrorCode DMDACreate(MPI_Comm comm, DM *da)
463: {
465:   DMCreate(comm, da);
466:   DMSetType(*da, DMDA);
467:   return 0;
468: }