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