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