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