xref: /petsc/src/dm/impls/da/dacreate.c (revision 1ebf93c6b7d760d592de6ebe6cdc0debaa3caf75)
1 
2 #include <petsc/private/dmdaimpl.h>    /*I   "petscdmda.h"   I*/
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "DMSetFromOptions_DA"
6 PetscErrorCode  DMSetFromOptions_DA(PetscOptionItems *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 = PetscOptionsIntArray("-da_refine_hierarchy_x","Refinement factor for each level","None",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 = PetscOptionsIntArray("-da_refine_hierarchy_y","Refinement factor for each level","None",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 = PetscOptionsIntArray("-da_refine_hierarchy_z","Refinement factor for each level","None",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 extern PetscErrorCode DMProjectFunctionLocal_DA(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
335 extern PetscErrorCode DMComputeL2Diff_DA(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
336 extern PetscErrorCode DMComputeL2GradientDiff_DA(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal [], const PetscReal [],PetscInt, PetscScalar *, void *), void **, Vec,const PetscReal [], PetscReal *);
337 
338 
339 #undef __FUNCT__
340 #define __FUNCT__ "DMCreate_DA"
341 PETSC_EXTERN PetscErrorCode DMCreate_DA(DM da)
342 {
343   PetscErrorCode ierr;
344   DM_DA          *dd;
345 
346   PetscFunctionBegin;
347   PetscValidPointer(da,1);
348   ierr     = PetscNewLog(da,&dd);CHKERRQ(ierr);
349   da->data = dd;
350 
351   da->dim        = -1;
352   dd->interptype = DMDA_Q1;
353   dd->refine_x   = 2;
354   dd->refine_y   = 2;
355   dd->refine_z   = 2;
356   dd->coarsen_x  = 2;
357   dd->coarsen_y  = 2;
358   dd->coarsen_z  = 2;
359   dd->fieldname  = NULL;
360   dd->nlocal     = -1;
361   dd->Nlocal     = -1;
362   dd->M          = -1;
363   dd->N          = -1;
364   dd->P          = -1;
365   dd->m          = -1;
366   dd->n          = -1;
367   dd->p          = -1;
368   dd->w          = -1;
369   dd->s          = -1;
370 
371   dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1;
372   dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1;
373 
374   dd->Nsub            = 1;
375   dd->xol             = 0;
376   dd->yol             = 0;
377   dd->zol             = 0;
378   dd->xo              = 0;
379   dd->yo              = 0;
380   dd->zo              = 0;
381   dd->Mo              = -1;
382   dd->No              = -1;
383   dd->Po              = -1;
384 
385   dd->gtol         = NULL;
386   dd->ltol         = NULL;
387   dd->ao           = NULL;
388   PetscStrallocpy(AOBASIC,(char**)&dd->aotype);
389   dd->base         = -1;
390   dd->bx           = DM_BOUNDARY_NONE;
391   dd->by           = DM_BOUNDARY_NONE;
392   dd->bz           = DM_BOUNDARY_NONE;
393   dd->stencil_type = DMDA_STENCIL_BOX;
394   dd->interptype   = DMDA_Q1;
395   dd->lx           = NULL;
396   dd->ly           = NULL;
397   dd->lz           = NULL;
398 
399   dd->elementtype = DMDA_ELEMENT_Q1;
400 
401   da->ops->globaltolocalbegin          = DMGlobalToLocalBegin_DA;
402   da->ops->globaltolocalend            = DMGlobalToLocalEnd_DA;
403   da->ops->localtoglobalbegin          = DMLocalToGlobalBegin_DA;
404   da->ops->localtoglobalend            = DMLocalToGlobalEnd_DA;
405   da->ops->localtolocalbegin           = DMLocalToLocalBegin_DA;
406   da->ops->localtolocalend             = DMLocalToLocalEnd_DA;
407   da->ops->createglobalvector          = DMCreateGlobalVector_DA;
408   da->ops->createlocalvector           = DMCreateLocalVector_DA;
409   da->ops->createinterpolation         = DMCreateInterpolation_DA;
410   da->ops->getcoloring                 = DMCreateColoring_DA;
411   da->ops->creatematrix                = DMCreateMatrix_DA;
412   da->ops->refine                      = DMRefine_DA;
413   da->ops->coarsen                     = DMCoarsen_DA;
414   da->ops->refinehierarchy             = DMRefineHierarchy_DA;
415   da->ops->coarsenhierarchy            = DMCoarsenHierarchy_DA;
416   da->ops->getinjection                = DMCreateInjection_DA;
417   da->ops->getaggregates               = DMCreateAggregates_DA;
418   da->ops->destroy                     = DMDestroy_DA;
419   da->ops->view                        = 0;
420   da->ops->setfromoptions              = DMSetFromOptions_DA;
421   da->ops->setup                       = DMSetUp_DA;
422   da->ops->clone                       = DMClone_DA;
423   da->ops->load                        = DMLoad_DA;
424   da->ops->createcoordinatedm          = DMCreateCoordinateDM_DA;
425   da->ops->createsubdm                 = DMCreateSubDM_DA;
426   da->ops->createfielddecomposition    = DMCreateFieldDecomposition_DA;
427   da->ops->createdomaindecomposition   = DMCreateDomainDecomposition_DA;
428   da->ops->createddscatters            = DMCreateDomainDecompositionScatters_DA;
429   da->ops->getdimpoints                = DMGetDimPoints_DA;
430   da->ops->projectfunctionlocal        = DMProjectFunctionLocal_DA;
431   da->ops->computel2diff               = DMComputeL2Diff_DA;
432   da->ops->computel2gradientdiff       = DMComputeL2GradientDiff_DA;
433   PetscFunctionReturn(0);
434 }
435 
436 #undef __FUNCT__
437 #define __FUNCT__ "DMDACreate"
438 /*@
439   DMDACreate - Creates a DMDA object.
440 
441   Collective on MPI_Comm
442 
443   Input Parameter:
444 . comm - The communicator for the DMDA object
445 
446   Output Parameter:
447 . da  - The DMDA object
448 
449   Level: advanced
450 
451   Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist?
452 
453 .keywords: DMDA, create
454 .seealso:  DMDASetSizes(), DMDADuplicate(),  DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
455 @*/
456 PetscErrorCode  DMDACreate(MPI_Comm comm, DM *da)
457 {
458   PetscErrorCode ierr;
459 
460   PetscFunctionBegin;
461   PetscValidPointer(da,2);
462   ierr = DMCreate(comm,da);CHKERRQ(ierr);
463   ierr = DMSetType(*da,DMDA);CHKERRQ(ierr);
464   PetscFunctionReturn(0);
465 }
466 
467 
468