xref: /petsc/src/dm/impls/da/dacreate.c (revision bebe2cf65d55febe21a5af8db2bd2e168caaa2e7)
1 
2 #include <petsc/private/dmdaimpl.h>    /*I   "petscdmda.h"   I*/
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "DMSetFromOptions_DA"
6 PetscErrorCode  DMSetFromOptions_DA(PetscOptions *PetscOptionsObject,DM da)
7 {
8   PetscErrorCode ierr;
9   DM_DA          *dd    = (DM_DA*)da->data;
10   PetscInt       refine = 0,dim = da->dim,maxnlevels = 100,refx[100],refy[100],refz[100],n,i;
11   PetscBool      bM     = PETSC_FALSE,bN = PETSC_FALSE, bP = PETSC_FALSE,flg;
12 
13   PetscFunctionBegin;
14   PetscValidHeaderSpecific(da,DM_CLASSID,1);
15 
16   if (dd->M < 0) {
17     dd->M           = -dd->M;
18     bM              = PETSC_TRUE;
19     dd->negativeMNP = PETSC_TRUE;
20   }
21   if (da->dim > 1 && dd->N < 0) {
22     dd->N           = -dd->N;
23     bN              = PETSC_TRUE;
24     dd->negativeMNP = PETSC_TRUE;
25   }
26   if (da->dim > 2 && dd->P < 0) {
27     dd->P           = -dd->P;
28     bP              = PETSC_TRUE;
29     dd->negativeMNP = PETSC_TRUE;
30   }
31 
32   ierr = PetscOptionsHead(PetscOptionsObject,"DMDA Options");CHKERRQ(ierr);
33   if (bM) {ierr = PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DMDASetSizes",dd->M,&dd->M,NULL);CHKERRQ(ierr);}
34   if (bN) {ierr = PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DMDASetSizes",dd->N,&dd->N,NULL);CHKERRQ(ierr);}
35   if (bP) {ierr = PetscOptionsInt("-da_grid_z","Number of grid points in z direction","DMDASetSizes",dd->P,&dd->P,NULL);CHKERRQ(ierr);}
36 
37   ierr = PetscOptionsInt("-da_overlap","Decomposition overlap in all directions","DMDASetOverlap",dd->xol,&dd->xol,&flg);CHKERRQ(ierr);
38   if (flg) {ierr = DMDASetOverlap(da,dd->xol,dd->xol,dd->xol);CHKERRQ(ierr);}
39   ierr = PetscOptionsInt("-da_overlap_x","Decomposition overlap in x direction","DMDASetOverlap",dd->xol,&dd->xol,NULL);CHKERRQ(ierr);
40   if (dim > 1) {ierr = PetscOptionsInt("-da_overlap_y","Decomposition overlap in y direction","DMDASetOverlap",dd->yol,&dd->yol,NULL);CHKERRQ(ierr);}
41   if (dim > 2) {ierr = PetscOptionsInt("-da_overlap_z","Decomposition overlap in z direction","DMDASetOverlap",dd->zol,&dd->zol,NULL);CHKERRQ(ierr);}
42 
43   ierr = PetscOptionsInt("-da_local_subdomains","","DMDASetNumLocalSubdomains",dd->Nsub,&dd->Nsub,&flg);CHKERRQ(ierr);
44   if (flg) {ierr = DMDASetNumLocalSubDomains(da,dd->Nsub);CHKERRQ(ierr);}
45 
46   /* Handle DMDA parallel distibution */
47   ierr = PetscOptionsInt("-da_processors_x","Number of processors in x direction","DMDASetNumProcs",dd->m,&dd->m,NULL);CHKERRQ(ierr);
48   if (dim > 1) {ierr = PetscOptionsInt("-da_processors_y","Number of processors in y direction","DMDASetNumProcs",dd->n,&dd->n,NULL);CHKERRQ(ierr);}
49   if (dim > 2) {ierr = PetscOptionsInt("-da_processors_z","Number of processors in z direction","DMDASetNumProcs",dd->p,&dd->p,NULL);CHKERRQ(ierr);}
50   /* Handle DMDA refinement */
51   ierr = PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DMDASetRefinementFactor",dd->refine_x,&dd->refine_x,NULL);CHKERRQ(ierr);
52   if (dim > 1) {ierr = PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DMDASetRefinementFactor",dd->refine_y,&dd->refine_y,NULL);CHKERRQ(ierr);}
53   if (dim > 2) {ierr = PetscOptionsInt("-da_refine_z","Refinement ratio in z direction","DMDASetRefinementFactor",dd->refine_z,&dd->refine_z,NULL);CHKERRQ(ierr);}
54   dd->coarsen_x = dd->refine_x; dd->coarsen_y = dd->refine_y; dd->coarsen_z = dd->refine_z;
55 
56   /* Get refinement factors, defaults taken from the coarse DMDA */
57   ierr = DMDAGetRefinementFactor(da,&refx[0],&refy[0],&refz[0]);CHKERRQ(ierr);
58   for (i=1; i<maxnlevels; i++) {
59     refx[i] = refx[0];
60     refy[i] = refy[0];
61     refz[i] = refz[0];
62   }
63   n    = maxnlevels;
64   ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);CHKERRQ(ierr);
65   if (da->levelup - da->leveldown >= 0) dd->refine_x = refx[da->levelup - da->leveldown];
66   if (da->levelup - da->leveldown >= 1) dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
67   if (dim > 1) {
68     n    = maxnlevels;
69     ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);CHKERRQ(ierr);
70     if (da->levelup - da->leveldown >= 0) dd->refine_y = refy[da->levelup - da->leveldown];
71     if (da->levelup - da->leveldown >= 1) dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
72   }
73   if (dim > 2) {
74     n    = maxnlevels;
75     ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);CHKERRQ(ierr);
76     if (da->levelup - da->leveldown >= 0) dd->refine_z = refz[da->levelup - da->leveldown];
77     if (da->levelup - da->leveldown >= 1) dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
78   }
79 
80   if (dd->negativeMNP) {ierr = PetscOptionsInt("-da_refine","Uniformly refine DA one or more times","None",refine,&refine,NULL);CHKERRQ(ierr);}
81   ierr = PetscOptionsTail();CHKERRQ(ierr);
82 
83   while (refine--) {
84     if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
85       dd->M = dd->refine_x*dd->M;
86     } else {
87       dd->M = 1 + dd->refine_x*(dd->M - 1);
88     }
89     if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
90       dd->N = dd->refine_y*dd->N;
91     } else {
92       dd->N = 1 + dd->refine_y*(dd->N - 1);
93     }
94     if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
95       dd->P = dd->refine_z*dd->P;
96     } else {
97       dd->P = 1 + dd->refine_z*(dd->P - 1);
98     }
99     da->levelup++;
100     if (da->levelup - da->leveldown >= 0) {
101       dd->refine_x = refx[da->levelup - da->leveldown];
102       dd->refine_y = refy[da->levelup - da->leveldown];
103       dd->refine_z = refz[da->levelup - da->leveldown];
104     }
105     if (da->levelup - da->leveldown >= 1) {
106       dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
107       dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
108       dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
109     }
110   }
111   PetscFunctionReturn(0);
112 }
113 
114 extern PetscErrorCode  DMCreateGlobalVector_DA(DM,Vec*);
115 extern PetscErrorCode  DMCreateLocalVector_DA(DM,Vec*);
116 extern PetscErrorCode  DMGlobalToLocalBegin_DA(DM,Vec,InsertMode,Vec);
117 extern PetscErrorCode  DMGlobalToLocalEnd_DA(DM,Vec,InsertMode,Vec);
118 extern PetscErrorCode  DMLocalToGlobalBegin_DA(DM,Vec,InsertMode,Vec);
119 extern PetscErrorCode  DMLocalToGlobalEnd_DA(DM,Vec,InsertMode,Vec);
120 extern PetscErrorCode  DMLocalToLocalBegin_DA(DM,Vec,InsertMode,Vec);
121 extern PetscErrorCode  DMLocalToLocalEnd_DA(DM,Vec,InsertMode,Vec);
122 extern PetscErrorCode  DMCreateInterpolation_DA(DM,DM,Mat*,Vec*);
123 extern PetscErrorCode  DMCreateColoring_DA(DM,ISColoringType,ISColoring*);
124 extern PetscErrorCode  DMCreateMatrix_DA(DM,Mat*);
125 extern PetscErrorCode  DMCreateCoordinateDM_DA(DM,DM*);
126 extern PetscErrorCode  DMRefine_DA(DM,MPI_Comm,DM*);
127 extern PetscErrorCode  DMCoarsen_DA(DM,MPI_Comm,DM*);
128 extern PetscErrorCode  DMRefineHierarchy_DA(DM,PetscInt,DM[]);
129 extern PetscErrorCode  DMCoarsenHierarchy_DA(DM,PetscInt,DM[]);
130 extern PetscErrorCode  DMCreateInjection_DA(DM,DM,Mat*);
131 extern PetscErrorCode  DMCreateAggregates_DA(DM,DM,Mat*);
132 extern PetscErrorCode  DMView_DA(DM,PetscViewer);
133 extern PetscErrorCode  DMSetUp_DA(DM);
134 extern PetscErrorCode  DMDestroy_DA(DM);
135 extern PetscErrorCode  DMCreateDomainDecomposition_DA(DM,PetscInt*,char***,IS**,IS**,DM**);
136 extern PetscErrorCode  DMCreateDomainDecompositionScatters_DA(DM,PetscInt,DM*,VecScatter**,VecScatter**,VecScatter**);
137 
138 #undef __FUNCT__
139 #define __FUNCT__ "DMLoad_DA"
140 PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer)
141 {
142   PetscErrorCode   ierr;
143   PetscInt         dim,m,n,p,dof,swidth;
144   DMDAStencilType  stencil;
145   DMBoundaryType   bx,by,bz;
146   PetscBool        coors;
147   DM               dac;
148   Vec              c;
149 
150   PetscFunctionBegin;
151   ierr = PetscViewerBinaryRead(viewer,&dim,1,NULL,PETSC_INT);CHKERRQ(ierr);
152   ierr = PetscViewerBinaryRead(viewer,&m,1,NULL,PETSC_INT);CHKERRQ(ierr);
153   ierr = PetscViewerBinaryRead(viewer,&n,1,NULL,PETSC_INT);CHKERRQ(ierr);
154   ierr = PetscViewerBinaryRead(viewer,&p,1,NULL,PETSC_INT);CHKERRQ(ierr);
155   ierr = PetscViewerBinaryRead(viewer,&dof,1,NULL,PETSC_INT);CHKERRQ(ierr);
156   ierr = PetscViewerBinaryRead(viewer,&swidth,1,NULL,PETSC_INT);CHKERRQ(ierr);
157   ierr = PetscViewerBinaryRead(viewer,&bx,1,NULL,PETSC_ENUM);CHKERRQ(ierr);
158   ierr = PetscViewerBinaryRead(viewer,&by,1,NULL,PETSC_ENUM);CHKERRQ(ierr);
159   ierr = PetscViewerBinaryRead(viewer,&bz,1,NULL,PETSC_ENUM);CHKERRQ(ierr);
160   ierr = PetscViewerBinaryRead(viewer,&stencil,1,NULL,PETSC_ENUM);CHKERRQ(ierr);
161 
162   ierr = DMSetDimension(da, dim);CHKERRQ(ierr);
163   ierr = DMDASetSizes(da, m,n,p);CHKERRQ(ierr);
164   ierr = DMDASetBoundaryType(da, bx, by, bz);CHKERRQ(ierr);
165   ierr = DMDASetDof(da, dof);CHKERRQ(ierr);
166   ierr = DMDASetStencilType(da, stencil);CHKERRQ(ierr);
167   ierr = DMDASetStencilWidth(da, swidth);CHKERRQ(ierr);
168   ierr = DMSetUp(da);CHKERRQ(ierr);
169   ierr = PetscViewerBinaryRead(viewer,&coors,1,NULL,PETSC_ENUM);CHKERRQ(ierr);
170   if (coors) {
171     ierr = DMGetCoordinateDM(da,&dac);CHKERRQ(ierr);
172     ierr = DMCreateGlobalVector(dac,&c);CHKERRQ(ierr);
173     ierr = VecLoad(c,viewer);CHKERRQ(ierr);
174     ierr = DMSetCoordinates(da,c);CHKERRQ(ierr);
175     ierr = VecDestroy(&c);CHKERRQ(ierr);
176   }
177   PetscFunctionReturn(0);
178 }
179 
180 #undef __FUNCT__
181 #define __FUNCT__ "DMCreateSubDM_DA"
182 PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
183 {
184   DM_DA         *da = (DM_DA*) dm->data;
185   PetscSection   section;
186   PetscErrorCode ierr;
187 
188   PetscFunctionBegin;
189   if (subdm) {
190     PetscSF sf;
191     Vec     coords;
192     void   *ctx;
193     /* Cannot use DMClone since the dof stuff is mixed in. Ugh
194     ierr = DMClone(dm, subdm);CHKERRQ(ierr); */
195     ierr = DMCreate(PetscObjectComm((PetscObject)dm), subdm);CHKERRQ(ierr);
196     ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
197     ierr = DMSetPointSF(*subdm, sf);CHKERRQ(ierr);
198     ierr = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr);
199     ierr = DMSetApplicationContext(*subdm, ctx);CHKERRQ(ierr);
200     ierr = DMGetCoordinatesLocal(dm, &coords);CHKERRQ(ierr);
201     if (coords) {
202       ierr = DMSetCoordinatesLocal(*subdm, coords);CHKERRQ(ierr);
203     } else {
204       ierr = DMGetCoordinates(dm, &coords);CHKERRQ(ierr);
205       if (coords) {ierr = DMSetCoordinates(*subdm, coords);CHKERRQ(ierr);}
206     }
207 
208     ierr = DMSetType(*subdm, DMDA);CHKERRQ(ierr);
209     ierr = DMSetDimension(*subdm, dm->dim);CHKERRQ(ierr);
210     ierr = DMDASetSizes(*subdm, da->M, da->N, da->P);CHKERRQ(ierr);
211     ierr = DMDASetNumProcs(*subdm, da->m, da->n, da->p);CHKERRQ(ierr);
212     ierr = DMDASetBoundaryType(*subdm, da->bx, da->by, da->bz);CHKERRQ(ierr);
213     ierr = DMDASetDof(*subdm, numFields);CHKERRQ(ierr);
214     ierr = DMDASetStencilType(*subdm, da->stencil_type);CHKERRQ(ierr);
215     ierr = DMDASetStencilWidth(*subdm, da->s);CHKERRQ(ierr);
216     ierr = DMDASetOwnershipRanges(*subdm, da->lx, da->ly, da->lz);CHKERRQ(ierr);
217   }
218   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
219   if (section) {
220     ierr = DMCreateSubDM_Section_Private(dm, numFields, fields, is, subdm);CHKERRQ(ierr);
221   } else {
222     if (is) {
223       PetscInt *indices, cnt = 0, dof = da->w, i, j;
224 
225       ierr = PetscMalloc1(da->Nlocal*numFields/dof, &indices);CHKERRQ(ierr);
226       for (i = da->base/dof; i < (da->base+da->Nlocal)/dof; ++i) {
227         for (j = 0; j < numFields; ++j) {
228           indices[cnt++] = dof*i + fields[j];
229         }
230       }
231       if (cnt != da->Nlocal*numFields/dof) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count %d does not equal expected value %d", cnt, da->Nlocal*numFields/dof);
232       ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), cnt, indices, PETSC_OWN_POINTER, is);CHKERRQ(ierr);
233     }
234   }
235   PetscFunctionReturn(0);
236 }
237 
238 #undef __FUNCT__
239 #define __FUNCT__ "DMCreateFieldDecomposition_DA"
240 PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
241 {
242   PetscInt       i;
243   PetscErrorCode ierr;
244   DM_DA          *dd = (DM_DA*)dm->data;
245   PetscInt       dof = dd->w;
246 
247   PetscFunctionBegin;
248   if (len) *len = dof;
249   if (islist) {
250     Vec      v;
251     PetscInt rstart,n;
252 
253     ierr = DMGetGlobalVector(dm,&v);CHKERRQ(ierr);
254     ierr = VecGetOwnershipRange(v,&rstart,NULL);CHKERRQ(ierr);
255     ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
256     ierr = DMRestoreGlobalVector(dm,&v);CHKERRQ(ierr);
257     ierr = PetscMalloc1(dof,islist);CHKERRQ(ierr);
258     for (i=0; i<dof; i++) {
259       ierr = ISCreateStride(PetscObjectComm((PetscObject)dm),n/dof,rstart+i,dof,&(*islist)[i]);CHKERRQ(ierr);
260     }
261   }
262   if (namelist) {
263     ierr = PetscMalloc1(dof, namelist);CHKERRQ(ierr);
264     if (dd->fieldname) {
265       for (i=0; i<dof; i++) {
266         ierr = PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);CHKERRQ(ierr);
267       }
268     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames");
269   }
270   if (dmlist) {
271     DM da;
272 
273     ierr = DMDACreate(PetscObjectComm((PetscObject)dm), &da);CHKERRQ(ierr);
274     ierr = DMSetDimension(da, dm->dim);CHKERRQ(ierr);
275     ierr = DMDASetSizes(da, dd->M, dd->N, dd->P);CHKERRQ(ierr);
276     ierr = DMDASetNumProcs(da, dd->m, dd->n, dd->p);CHKERRQ(ierr);
277     ierr = DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);CHKERRQ(ierr);
278     ierr = DMDASetDof(da, 1);CHKERRQ(ierr);
279     ierr = DMDASetStencilType(da, dd->stencil_type);CHKERRQ(ierr);
280     ierr = DMDASetStencilWidth(da, dd->s);CHKERRQ(ierr);
281     ierr = DMSetUp(da);CHKERRQ(ierr);
282     ierr = PetscMalloc1(dof,dmlist);CHKERRQ(ierr);
283     for (i=0; i<dof-1; i++) {ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);}
284     for (i=0; i<dof; i++) (*dmlist)[i] = da;
285   }
286   PetscFunctionReturn(0);
287 }
288 
289 #undef __FUNCT__
290 #define __FUNCT__ "DMClone_DA"
291 PetscErrorCode DMClone_DA(DM dm, DM *newdm)
292 {
293   DM_DA         *da = (DM_DA *) dm->data;
294   PetscErrorCode ierr;
295 
296   PetscFunctionBegin;
297   ierr = DMSetType(*newdm, DMDA);CHKERRQ(ierr);
298   ierr = DMSetDimension(*newdm, dm->dim);CHKERRQ(ierr);
299   ierr = DMDASetSizes(*newdm, da->M, da->N, da->P);CHKERRQ(ierr);
300   ierr = DMDASetNumProcs(*newdm, da->m, da->n, da->p);CHKERRQ(ierr);
301   ierr = DMDASetBoundaryType(*newdm, da->bx, da->by, da->bz);CHKERRQ(ierr);
302   ierr = DMDASetDof(*newdm, da->w);CHKERRQ(ierr);
303   ierr = DMDASetStencilType(*newdm, da->stencil_type);CHKERRQ(ierr);
304   ierr = DMDASetStencilWidth(*newdm, da->s);CHKERRQ(ierr);
305   ierr = DMDASetOwnershipRanges(*newdm, da->lx, da->ly, da->lz);CHKERRQ(ierr);
306   ierr = DMSetUp(*newdm);CHKERRQ(ierr);
307   PetscFunctionReturn(0);
308 }
309 
310 #undef __FUNCT__
311 #define __FUNCT__ "DMGetDimPoints_DA"
312 static PetscErrorCode DMGetDimPoints_DA(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
313 {
314   PetscErrorCode ierr;
315 
316   PetscFunctionBegin;
317   ierr = DMDAGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);
318   PetscFunctionReturn(0);
319 }
320 
321 /*MC
322    DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions.
323          In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points.
324          In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width.
325 
326          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
327          vertex centered.
328 
329   Level: intermediate
330 
331 .seealso: DMType, DMCOMPOSITE, DMDACreate(), DMCreate(), DMSetType()
332 M*/
333 
334 
335 #undef __FUNCT__
336 #define __FUNCT__ "DMCreate_DA"
337 PETSC_EXTERN PetscErrorCode DMCreate_DA(DM da)
338 {
339   PetscErrorCode ierr;
340   DM_DA          *dd;
341 
342   PetscFunctionBegin;
343   PetscValidPointer(da,1);
344   ierr     = PetscNewLog(da,&dd);CHKERRQ(ierr);
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; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1;
368   dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1;
369 
370   dd->Nsub            = 1;
371   dd->xol             = 0;
372   dd->yol             = 0;
373   dd->zol             = 0;
374   dd->xo              = 0;
375   dd->yo              = 0;
376   dd->zo              = 0;
377   dd->Mo              = -1;
378   dd->No              = -1;
379   dd->Po              = -1;
380 
381   dd->gtol         = NULL;
382   dd->ltol         = NULL;
383   dd->ao           = NULL;
384   PetscStrallocpy(AOBASIC,(char**)&dd->aotype);
385   dd->base         = -1;
386   dd->bx           = DM_BOUNDARY_NONE;
387   dd->by           = DM_BOUNDARY_NONE;
388   dd->bz           = DM_BOUNDARY_NONE;
389   dd->stencil_type = DMDA_STENCIL_BOX;
390   dd->interptype   = DMDA_Q1;
391   dd->lx           = NULL;
392   dd->ly           = NULL;
393   dd->lz           = NULL;
394 
395   dd->elementtype = DMDA_ELEMENT_Q1;
396 
397   da->ops->globaltolocalbegin          = DMGlobalToLocalBegin_DA;
398   da->ops->globaltolocalend            = DMGlobalToLocalEnd_DA;
399   da->ops->localtoglobalbegin          = DMLocalToGlobalBegin_DA;
400   da->ops->localtoglobalend            = DMLocalToGlobalEnd_DA;
401   da->ops->localtolocalbegin           = DMLocalToLocalBegin_DA;
402   da->ops->localtolocalend             = DMLocalToLocalEnd_DA;
403   da->ops->createglobalvector          = DMCreateGlobalVector_DA;
404   da->ops->createlocalvector           = DMCreateLocalVector_DA;
405   da->ops->createinterpolation         = DMCreateInterpolation_DA;
406   da->ops->getcoloring                 = DMCreateColoring_DA;
407   da->ops->creatematrix                = DMCreateMatrix_DA;
408   da->ops->refine                      = DMRefine_DA;
409   da->ops->coarsen                     = DMCoarsen_DA;
410   da->ops->refinehierarchy             = DMRefineHierarchy_DA;
411   da->ops->coarsenhierarchy            = DMCoarsenHierarchy_DA;
412   da->ops->getinjection                = DMCreateInjection_DA;
413   da->ops->getaggregates               = DMCreateAggregates_DA;
414   da->ops->destroy                     = DMDestroy_DA;
415   da->ops->view                        = 0;
416   da->ops->setfromoptions              = DMSetFromOptions_DA;
417   da->ops->setup                       = DMSetUp_DA;
418   da->ops->clone                       = DMClone_DA;
419   da->ops->load                        = DMLoad_DA;
420   da->ops->createcoordinatedm          = DMCreateCoordinateDM_DA;
421   da->ops->createsubdm                 = DMCreateSubDM_DA;
422   da->ops->createfielddecomposition    = DMCreateFieldDecomposition_DA;
423   da->ops->createdomaindecomposition   = DMCreateDomainDecomposition_DA;
424   da->ops->createddscatters            = DMCreateDomainDecompositionScatters_DA;
425   da->ops->getdimpoints                = DMGetDimPoints_DA;
426   PetscFunctionReturn(0);
427 }
428 
429 #undef __FUNCT__
430 #define __FUNCT__ "DMDACreate"
431 /*@
432   DMDACreate - Creates a DMDA object.
433 
434   Collective on MPI_Comm
435 
436   Input Parameter:
437 . comm - The communicator for the DMDA object
438 
439   Output Parameter:
440 . da  - The DMDA object
441 
442   Level: advanced
443 
444   Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist?
445 
446 .keywords: DMDA, create
447 .seealso:  DMDASetSizes(), DMDADuplicate(),  DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
448 @*/
449 PetscErrorCode  DMDACreate(MPI_Comm comm, DM *da)
450 {
451   PetscErrorCode ierr;
452 
453   PetscFunctionBegin;
454   PetscValidPointer(da,2);
455   ierr = DMCreate(comm,da);CHKERRQ(ierr);
456   ierr = DMSetType(*da,DMDA);CHKERRQ(ierr);
457   PetscFunctionReturn(0);
458 }
459 
460 
461