xref: /petsc/src/dm/impls/da/dacreate.c (revision a69119a591a03a9d906b29c0a4e9802e4d7c9795)
1 
2 #include <petsc/private/dmdaimpl.h> /*I   "petscdmda.h"   I*/
3 
4 PetscErrorCode DMSetFromOptions_DA(DM da, PetscOptionItems *PetscOptionsObject) {
5   DM_DA    *dd     = (DM_DA *)da->data;
6   PetscInt  refine = 0, dim = da->dim, maxnlevels = 100, refx[100], refy[100], refz[100], n, i;
7   PetscBool flg;
8 
9   PetscFunctionBegin;
10   PetscCheck(dd->M >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Dimension must be non-negative, call DMSetFromOptions() if you want to change the value at runtime");
11   PetscCheck(dd->N >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Dimension must be non-negative, call DMSetFromOptions() if you want to change the value at runtime");
12   PetscCheck(dd->P >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Dimension must be non-negative, call DMSetFromOptions() if you want to change the value at runtime");
13 
14   PetscOptionsHeadBegin(PetscOptionsObject, "DMDA Options");
15   PetscCall(PetscOptionsBoundedInt("-da_grid_x", "Number of grid points in x direction", "DMDASetSizes", dd->M, &dd->M, NULL, 1));
16   if (dim > 1) PetscCall(PetscOptionsBoundedInt("-da_grid_y", "Number of grid points in y direction", "DMDASetSizes", dd->N, &dd->N, NULL, 1));
17   if (dim > 2) PetscCall(PetscOptionsBoundedInt("-da_grid_z", "Number of grid points in z direction", "DMDASetSizes", dd->P, &dd->P, NULL, 1));
18 
19   PetscCall(PetscOptionsBoundedInt("-da_overlap", "Decomposition overlap in all directions", "DMDASetOverlap", dd->xol, &dd->xol, &flg, 0));
20   if (flg) PetscCall(DMDASetOverlap(da, dd->xol, dd->xol, dd->xol));
21   PetscCall(PetscOptionsBoundedInt("-da_overlap_x", "Decomposition overlap in x direction", "DMDASetOverlap", dd->xol, &dd->xol, NULL, 0));
22   if (dim > 1) PetscCall(PetscOptionsBoundedInt("-da_overlap_y", "Decomposition overlap in y direction", "DMDASetOverlap", dd->yol, &dd->yol, NULL, 0));
23   if (dim > 2) PetscCall(PetscOptionsBoundedInt("-da_overlap_z", "Decomposition overlap in z direction", "DMDASetOverlap", dd->zol, &dd->zol, NULL, 0));
24 
25   PetscCall(PetscOptionsBoundedInt("-da_local_subdomains", "", "DMDASetNumLocalSubdomains", dd->Nsub, &dd->Nsub, &flg, PETSC_DECIDE));
26   if (flg) PetscCall(DMDASetNumLocalSubDomains(da, dd->Nsub));
27 
28   /* Handle DMDA parallel distribution */
29   PetscCall(PetscOptionsBoundedInt("-da_processors_x", "Number of processors in x direction", "DMDASetNumProcs", dd->m, &dd->m, NULL, PETSC_DECIDE));
30   if (dim > 1) PetscCall(PetscOptionsBoundedInt("-da_processors_y", "Number of processors in y direction", "DMDASetNumProcs", dd->n, &dd->n, NULL, PETSC_DECIDE));
31   if (dim > 2) PetscCall(PetscOptionsBoundedInt("-da_processors_z", "Number of processors in z direction", "DMDASetNumProcs", dd->p, &dd->p, NULL, PETSC_DECIDE));
32   /* Handle DMDA refinement */
33   PetscCall(PetscOptionsBoundedInt("-da_refine_x", "Refinement ratio in x direction", "DMDASetRefinementFactor", dd->refine_x, &dd->refine_x, NULL, 1));
34   if (dim > 1) PetscCall(PetscOptionsBoundedInt("-da_refine_y", "Refinement ratio in y direction", "DMDASetRefinementFactor", dd->refine_y, &dd->refine_y, NULL, 1));
35   if (dim > 2) PetscCall(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;
39 
40   /* Get refinement factors, defaults taken from the coarse DMDA */
41   PetscCall(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   PetscCall(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     PetscCall(PetscMalloc1(n, &dd->refine_x_hier));
53     PetscCall(PetscArraycpy(dd->refine_x_hier, refx, n));
54   }
55   if (dim > 1) {
56     n = maxnlevels;
57     PetscCall(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       PetscCall(PetscMalloc1(n, &dd->refine_y_hier));
62       PetscCall(PetscArraycpy(dd->refine_y_hier, refy, n));
63     }
64   }
65   if (dim > 2) {
66     n = maxnlevels;
67     PetscCall(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       PetscCall(PetscMalloc1(n, &dd->refine_z_hier));
72       PetscCall(PetscArraycpy(dd->refine_z_hier, refz, n));
73     }
74   }
75 
76   PetscCall(PetscOptionsBoundedInt("-da_refine", "Uniformly refine DA one or more times", "None", refine, &refine, NULL, 0));
77   PetscOptionsHeadEnd();
78 
79   while (refine--) {
80     if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
81       PetscCall(PetscIntMultError(dd->refine_x, dd->M, &dd->M));
82     } else {
83       PetscCall(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       PetscCall(PetscIntMultError(dd->refine_y, dd->N, &dd->N));
88     } else {
89       PetscCall(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       PetscCall(PetscIntMultError(dd->refine_z, dd->P, &dd->P));
94     } else {
95       PetscCall(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   PetscFunctionReturn(0);
111 }
112 
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 *);
136 
137 PetscErrorCode DMLoad_DA(DM da, PetscViewer viewer) {
138   PetscInt        dim, m, n, p, dof, swidth;
139   DMDAStencilType stencil;
140   DMBoundaryType  bx, by, bz;
141   PetscBool       coors;
142   DM              dac;
143   Vec             c;
144 
145   PetscFunctionBegin;
146   PetscCall(PetscViewerBinaryRead(viewer, &dim, 1, NULL, PETSC_INT));
147   PetscCall(PetscViewerBinaryRead(viewer, &m, 1, NULL, PETSC_INT));
148   PetscCall(PetscViewerBinaryRead(viewer, &n, 1, NULL, PETSC_INT));
149   PetscCall(PetscViewerBinaryRead(viewer, &p, 1, NULL, PETSC_INT));
150   PetscCall(PetscViewerBinaryRead(viewer, &dof, 1, NULL, PETSC_INT));
151   PetscCall(PetscViewerBinaryRead(viewer, &swidth, 1, NULL, PETSC_INT));
152   PetscCall(PetscViewerBinaryRead(viewer, &bx, 1, NULL, PETSC_ENUM));
153   PetscCall(PetscViewerBinaryRead(viewer, &by, 1, NULL, PETSC_ENUM));
154   PetscCall(PetscViewerBinaryRead(viewer, &bz, 1, NULL, PETSC_ENUM));
155   PetscCall(PetscViewerBinaryRead(viewer, &stencil, 1, NULL, PETSC_ENUM));
156 
157   PetscCall(DMSetDimension(da, dim));
158   PetscCall(DMDASetSizes(da, m, n, p));
159   PetscCall(DMDASetBoundaryType(da, bx, by, bz));
160   PetscCall(DMDASetDof(da, dof));
161   PetscCall(DMDASetStencilType(da, stencil));
162   PetscCall(DMDASetStencilWidth(da, swidth));
163   PetscCall(DMSetUp(da));
164   PetscCall(PetscViewerBinaryRead(viewer, &coors, 1, NULL, PETSC_ENUM));
165   if (coors) {
166     PetscCall(DMGetCoordinateDM(da, &dac));
167     PetscCall(DMCreateGlobalVector(dac, &c));
168     PetscCall(VecLoad(c, viewer));
169     PetscCall(DMSetCoordinates(da, c));
170     PetscCall(VecDestroy(&c));
171   }
172   PetscFunctionReturn(0);
173 }
174 
175 PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm) {
176   DM_DA *da = (DM_DA *)dm->data;
177 
178   PetscFunctionBegin;
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     PetscCall(DMClone(dm, subdm)); */
185     PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm));
186     PetscCall(DMGetPointSF(dm, &sf));
187     PetscCall(DMSetPointSF(*subdm, sf));
188     PetscCall(DMGetApplicationContext(dm, &ctx));
189     PetscCall(DMSetApplicationContext(*subdm, ctx));
190     PetscCall(DMGetCoordinatesLocal(dm, &coords));
191     if (coords) {
192       PetscCall(DMSetCoordinatesLocal(*subdm, coords));
193     } else {
194       PetscCall(DMGetCoordinates(dm, &coords));
195       if (coords) PetscCall(DMSetCoordinates(*subdm, coords));
196     }
197 
198     PetscCall(DMSetType(*subdm, DMDA));
199     PetscCall(DMSetDimension(*subdm, dm->dim));
200     PetscCall(DMDASetSizes(*subdm, da->M, da->N, da->P));
201     PetscCall(DMDASetNumProcs(*subdm, da->m, da->n, da->p));
202     PetscCall(DMDASetBoundaryType(*subdm, da->bx, da->by, da->bz));
203     PetscCall(DMDASetDof(*subdm, numFields));
204     PetscCall(DMDASetStencilType(*subdm, da->stencil_type));
205     PetscCall(DMDASetStencilWidth(*subdm, da->s));
206     PetscCall(DMDASetOwnershipRanges(*subdm, da->lx, da->ly, da->lz));
207   }
208   if (is) {
209     PetscInt *indices, cnt = 0, dof = da->w, i, j;
210 
211     PetscCall(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     }
215     PetscCheck(cnt == da->Nlocal * numFields / dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count %" PetscInt_FMT " does not equal expected value %" PetscInt_FMT, cnt, da->Nlocal * numFields / dof);
216     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), cnt, indices, PETSC_OWN_POINTER, is));
217   }
218   PetscFunctionReturn(0);
219 }
220 
221 PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist) {
222   PetscInt i;
223   DM_DA   *dd  = (DM_DA *)dm->data;
224   PetscInt dof = dd->w;
225 
226   PetscFunctionBegin;
227   if (len) *len = dof;
228   if (islist) {
229     Vec      v;
230     PetscInt rstart, n;
231 
232     PetscCall(DMGetGlobalVector(dm, &v));
233     PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
234     PetscCall(VecGetLocalSize(v, &n));
235     PetscCall(DMRestoreGlobalVector(dm, &v));
236     PetscCall(PetscMalloc1(dof, islist));
237     for (i = 0; i < dof; i++) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), n / dof, rstart + i, dof, &(*islist)[i]));
238   }
239   if (namelist) {
240     PetscCall(PetscMalloc1(dof, namelist));
241     if (dd->fieldname) {
242       for (i = 0; i < dof; i++) PetscCall(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;
247 
248     PetscCall(DMDACreate(PetscObjectComm((PetscObject)dm), &da));
249     PetscCall(DMSetDimension(da, dm->dim));
250     PetscCall(DMDASetSizes(da, dd->M, dd->N, dd->P));
251     PetscCall(DMDASetNumProcs(da, dd->m, dd->n, dd->p));
252     PetscCall(DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz));
253     PetscCall(DMDASetDof(da, 1));
254     PetscCall(DMDASetStencilType(da, dd->stencil_type));
255     PetscCall(DMDASetStencilWidth(da, dd->s));
256     PetscCall(DMSetUp(da));
257     PetscCall(PetscMalloc1(dof, dmlist));
258     for (i = 0; i < dof - 1; i++) PetscCall(PetscObjectReference((PetscObject)da));
259     for (i = 0; i < dof; i++) (*dmlist)[i] = da;
260   }
261   PetscFunctionReturn(0);
262 }
263 
264 PetscErrorCode DMClone_DA(DM dm, DM *newdm) {
265   DM_DA *da = (DM_DA *)dm->data;
266 
267   PetscFunctionBegin;
268   PetscCall(DMSetType(*newdm, DMDA));
269   PetscCall(DMSetDimension(*newdm, dm->dim));
270   PetscCall(DMDASetSizes(*newdm, da->M, da->N, da->P));
271   PetscCall(DMDASetNumProcs(*newdm, da->m, da->n, da->p));
272   PetscCall(DMDASetBoundaryType(*newdm, da->bx, da->by, da->bz));
273   PetscCall(DMDASetDof(*newdm, da->w));
274   PetscCall(DMDASetStencilType(*newdm, da->stencil_type));
275   PetscCall(DMDASetStencilWidth(*newdm, da->s));
276   PetscCall(DMDASetOwnershipRanges(*newdm, da->lx, da->ly, da->lz));
277   PetscCall(DMSetUp(*newdm));
278   PetscFunctionReturn(0);
279 }
280 
281 static PetscErrorCode DMHasCreateInjection_DA(DM dm, PetscBool *flg) {
282   DM_DA *da = (DM_DA *)dm->data;
283 
284   PetscFunctionBegin;
285   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
286   PetscValidBoolPointer(flg, 2);
287   *flg = da->interptype == DMDA_Q1 ? PETSC_TRUE : PETSC_FALSE;
288   PetscFunctionReturn(0);
289 }
290 
291 static PetscErrorCode DMGetDimPoints_DA(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd) {
292   PetscFunctionBegin;
293   PetscCall(DMDAGetDepthStratum(dm, dim, pStart, pEnd));
294   PetscFunctionReturn(0);
295 }
296 
297 static PetscErrorCode DMGetNeighbors_DA(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[]) {
298   PetscInt        dim;
299   DMDAStencilType st;
300 
301   PetscFunctionBegin;
302   PetscCall(DMDAGetNeighbors(dm, ranks));
303   PetscCall(DMDAGetInfo(dm, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &st));
304 
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: break;
319   }
320   PetscFunctionReturn(0);
321 }
322 
323 /*MC
324    DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions.
325          In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points.
326          In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width.
327 
328          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
329          vertex centered; see the documentation for DMSTAG, a similar DM implementation which supports these staggered grids.
330 
331   Level: intermediate
332 
333 .seealso: `DMType`, `DMCOMPOSITE`, `DMSTAG`, `DMDACreate()`, `DMCreate()`, `DMSetType()`
334 M*/
335 
336 extern PetscErrorCode       DMLocatePoints_DA_Regular(DM, Vec, DMPointLocationType, PetscSF);
337 PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_DMDA(PetscObject, PetscViewer);
338 
339 PETSC_EXTERN PetscErrorCode DMCreate_DA(DM da) {
340   DM_DA *dd;
341 
342   PetscFunctionBegin;
343   PetscValidPointer(da, 1);
344   PetscCall(PetscNewLog(da, &dd));
345   da->data = dd;
346 
347   da->dim        = -1;
348   dd->interptype = DMDA_Q1;
349   dd->refine_x   = 2;
350   dd->refine_y   = 2;
351   dd->refine_z   = 2;
352   dd->coarsen_x  = 2;
353   dd->coarsen_y  = 2;
354   dd->coarsen_z  = 2;
355   dd->fieldname  = NULL;
356   dd->nlocal     = -1;
357   dd->Nlocal     = -1;
358   dd->M          = -1;
359   dd->N          = -1;
360   dd->P          = -1;
361   dd->m          = -1;
362   dd->n          = -1;
363   dd->p          = -1;
364   dd->w          = -1;
365   dd->s          = -1;
366 
367   dd->xs = -1;
368   dd->xe = -1;
369   dd->ys = -1;
370   dd->ye = -1;
371   dd->zs = -1;
372   dd->ze = -1;
373   dd->Xs = -1;
374   dd->Xe = -1;
375   dd->Ys = -1;
376   dd->Ye = -1;
377   dd->Zs = -1;
378   dd->Ze = -1;
379 
380   dd->Nsub = 1;
381   dd->xol  = 0;
382   dd->yol  = 0;
383   dd->zol  = 0;
384   dd->xo   = 0;
385   dd->yo   = 0;
386   dd->zo   = 0;
387   dd->Mo   = -1;
388   dd->No   = -1;
389   dd->Po   = -1;
390 
391   dd->gtol = NULL;
392   dd->ltol = NULL;
393   dd->ao   = NULL;
394   PetscStrallocpy(AOBASIC, (char **)&dd->aotype);
395   dd->base         = -1;
396   dd->bx           = DM_BOUNDARY_NONE;
397   dd->by           = DM_BOUNDARY_NONE;
398   dd->bz           = DM_BOUNDARY_NONE;
399   dd->stencil_type = DMDA_STENCIL_BOX;
400   dd->interptype   = DMDA_Q1;
401   dd->lx           = NULL;
402   dd->ly           = NULL;
403   dd->lz           = NULL;
404 
405   dd->elementtype = DMDA_ELEMENT_Q1;
406 
407   da->ops->globaltolocalbegin        = DMGlobalToLocalBegin_DA;
408   da->ops->globaltolocalend          = DMGlobalToLocalEnd_DA;
409   da->ops->localtoglobalbegin        = DMLocalToGlobalBegin_DA;
410   da->ops->localtoglobalend          = DMLocalToGlobalEnd_DA;
411   da->ops->localtolocalbegin         = DMLocalToLocalBegin_DA;
412   da->ops->localtolocalend           = DMLocalToLocalEnd_DA;
413   da->ops->createglobalvector        = DMCreateGlobalVector_DA;
414   da->ops->createlocalvector         = DMCreateLocalVector_DA;
415   da->ops->createinterpolation       = DMCreateInterpolation_DA;
416   da->ops->getcoloring               = DMCreateColoring_DA;
417   da->ops->creatematrix              = DMCreateMatrix_DA;
418   da->ops->refine                    = DMRefine_DA;
419   da->ops->coarsen                   = DMCoarsen_DA;
420   da->ops->refinehierarchy           = DMRefineHierarchy_DA;
421   da->ops->coarsenhierarchy          = DMCoarsenHierarchy_DA;
422   da->ops->createinjection           = DMCreateInjection_DA;
423   da->ops->hascreateinjection        = DMHasCreateInjection_DA;
424   da->ops->destroy                   = DMDestroy_DA;
425   da->ops->view                      = NULL;
426   da->ops->setfromoptions            = DMSetFromOptions_DA;
427   da->ops->setup                     = DMSetUp_DA;
428   da->ops->clone                     = DMClone_DA;
429   da->ops->load                      = DMLoad_DA;
430   da->ops->createcoordinatedm        = DMCreateCoordinateDM_DA;
431   da->ops->createsubdm               = DMCreateSubDM_DA;
432   da->ops->createfielddecomposition  = DMCreateFieldDecomposition_DA;
433   da->ops->createdomaindecomposition = DMCreateDomainDecomposition_DA;
434   da->ops->createddscatters          = DMCreateDomainDecompositionScatters_DA;
435   da->ops->getdimpoints              = DMGetDimPoints_DA;
436   da->ops->getneighbors              = DMGetNeighbors_DA;
437   da->ops->locatepoints              = DMLocatePoints_DA_Regular;
438   da->ops->getcompatibility          = DMGetCompatibility_DA;
439   PetscCall(PetscObjectComposeFunction((PetscObject)da, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_DMDA));
440   PetscFunctionReturn(0);
441 }
442 
443 /*@
444   DMDACreate - Creates a DMDA object.
445 
446   Collective
447 
448   Input Parameter:
449 . comm - The communicator for the DMDA object
450 
451   Output Parameter:
452 . da  - The DMDA object
453 
454   Level: advanced
455 
456   Developers Note:
457   Since there exists DMDACreate1/2/3d() should this routine even exist?
458 
459 .seealso: `DMDASetSizes()`, `DMClone()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
460 @*/
461 PetscErrorCode DMDACreate(MPI_Comm comm, DM *da) {
462   PetscFunctionBegin;
463   PetscValidPointer(da, 2);
464   PetscCall(DMCreate(comm, da));
465   PetscCall(DMSetType(*da, DMDA));
466   PetscFunctionReturn(0);
467 }
468