xref: /petsc/src/dm/impls/da/dalocal.c (revision 3385ff7e2d1212a0f59ed2e0b7938fe2031de2c1) 
1 
2 /*
3   Code for manipulating distributed regular arrays in parallel.
4 */
5 
6 #include <petsc-private/dmdaimpl.h>    /*I   "petscdmda.h"   I*/
7 #include <petscsf.h>
8 
9 /*
10    This allows the DMDA vectors to properly tell MATLAB their dimensions
11 */
12 #if defined(PETSC_HAVE_MATLAB_ENGINE)
13 #include <engine.h>   /* MATLAB include file */
14 #include <mex.h>      /* MATLAB include file */
15 #undef __FUNCT__
16 #define __FUNCT__ "VecMatlabEnginePut_DA2d"
17 static PetscErrorCode  VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
18 {
19   PetscErrorCode ierr;
20   PetscInt       n,m;
21   Vec            vec = (Vec)obj;
22   PetscScalar    *array;
23   mxArray        *mat;
24   DM             da;
25 
26   PetscFunctionBegin;
27   ierr = VecGetDM(vec, &da);CHKERRQ(ierr);
28   if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA");
29   ierr = DMDAGetGhostCorners(da,0,0,0,&m,&n,0);CHKERRQ(ierr);
30 
31   ierr = VecGetArray(vec,&array);CHKERRQ(ierr);
32 #if !defined(PETSC_USE_COMPLEX)
33   mat = mxCreateDoubleMatrix(m,n,mxREAL);
34 #else
35   mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
36 #endif
37   ierr = PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));CHKERRQ(ierr);
38   ierr = PetscObjectName(obj);CHKERRQ(ierr);
39   engPutVariable((Engine*)mengine,obj->name,mat);
40 
41   ierr = VecRestoreArray(vec,&array);CHKERRQ(ierr);
42   PetscFunctionReturn(0);
43 }
44 #endif
45 
46 
47 #undef __FUNCT__
48 #define __FUNCT__ "DMCreateLocalVector_DA"
49 PetscErrorCode  DMCreateLocalVector_DA(DM da,Vec *g)
50 {
51   PetscErrorCode ierr;
52   DM_DA          *dd = (DM_DA*)da->data;
53 
54   PetscFunctionBegin;
55   PetscValidHeaderSpecific(da,DM_CLASSID,1);
56   PetscValidPointer(g,2);
57   if (da->defaultSection) {
58     ierr = DMCreateLocalVector_Section_Private(da,g);CHKERRQ(ierr);
59   } else {
60     ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr);
61     ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr);
62     ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr);
63     ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr);
64     ierr = VecSetDM(*g, da);CHKERRQ(ierr);
65 #if defined(PETSC_HAVE_MATLAB_ENGINE)
66     if (dd->w == 1  && dd->dim == 2) {
67       ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr);
68     }
69 #endif
70   }
71   PetscFunctionReturn(0);
72 }
73 
74 #undef __FUNCT__
75 #define __FUNCT__ "DMDAGetNumCells"
76 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
77 {
78   DM_DA          *da = (DM_DA*) dm->data;
79   const PetscInt dim = da->dim;
80   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
81   const PetscInt nC  = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
82 
83   PetscFunctionBegin;
84   if (numCellsX) {
85     PetscValidIntPointer(numCellsX,2);
86     *numCellsX = mx;
87   }
88   if (numCellsY) {
89     PetscValidIntPointer(numCellsX,3);
90     *numCellsY = my;
91   }
92   if (numCellsZ) {
93     PetscValidIntPointer(numCellsX,4);
94     *numCellsZ = mz;
95   }
96   if (numCells) {
97     PetscValidIntPointer(numCells,5);
98     *numCells = nC;
99   }
100   PetscFunctionReturn(0);
101 }
102 
103 #undef __FUNCT__
104 #define __FUNCT__ "DMDAGetNumVertices"
105 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
106 {
107   DM_DA          *da = (DM_DA*) dm->data;
108   const PetscInt dim = da->dim;
109   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
110   const PetscInt nVx = mx+1;
111   const PetscInt nVy = dim > 1 ? (my+1) : 1;
112   const PetscInt nVz = dim > 2 ? (mz+1) : 1;
113   const PetscInt nV  = nVx*nVy*nVz;
114 
115   PetscFunctionBegin;
116   if (numVerticesX) {
117     PetscValidIntPointer(numVerticesX,2);
118     *numVerticesX = nVx;
119   }
120   if (numVerticesY) {
121     PetscValidIntPointer(numVerticesY,3);
122     *numVerticesY = nVy;
123   }
124   if (numVerticesZ) {
125     PetscValidIntPointer(numVerticesZ,4);
126     *numVerticesZ = nVz;
127   }
128   if (numVertices) {
129     PetscValidIntPointer(numVertices,5);
130     *numVertices = nV;
131   }
132   PetscFunctionReturn(0);
133 }
134 
135 #undef __FUNCT__
136 #define __FUNCT__ "DMDAGetNumFaces"
137 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
138 {
139   DM_DA          *da = (DM_DA*) dm->data;
140   const PetscInt dim = da->dim;
141   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
142   const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
143   const PetscInt nXF = (mx+1)*nxF;
144   const PetscInt nyF = mx*(dim > 2 ? mz : 1);
145   const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0;
146   const PetscInt nzF = mx*(dim > 1 ? my : 0);
147   const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0;
148 
149   PetscFunctionBegin;
150   if (numXFacesX) {
151     PetscValidIntPointer(numXFacesX,2);
152     *numXFacesX = nxF;
153   }
154   if (numXFaces) {
155     PetscValidIntPointer(numXFaces,3);
156     *numXFaces = nXF;
157   }
158   if (numYFacesY) {
159     PetscValidIntPointer(numYFacesY,4);
160     *numYFacesY = nyF;
161   }
162   if (numYFaces) {
163     PetscValidIntPointer(numYFaces,5);
164     *numYFaces = nYF;
165   }
166   if (numZFacesZ) {
167     PetscValidIntPointer(numZFacesZ,6);
168     *numZFacesZ = nzF;
169   }
170   if (numZFaces) {
171     PetscValidIntPointer(numZFaces,7);
172     *numZFaces = nZF;
173   }
174   PetscFunctionReturn(0);
175 }
176 
177 #undef __FUNCT__
178 #define __FUNCT__ "DMDAGetHeightStratum"
179 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
180 {
181   DM_DA          *da = (DM_DA*) dm->data;
182   const PetscInt dim = da->dim;
183   PetscInt       nC, nV, nXF, nYF, nZF;
184   PetscErrorCode ierr;
185 
186   PetscFunctionBegin;
187   if (pStart) PetscValidIntPointer(pStart,3);
188   if (pEnd)   PetscValidIntPointer(pEnd,4);
189   ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
190   ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr);
191   ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr);
192   if (height == 0) {
193     /* Cells */
194     if (pStart) *pStart = 0;
195     if (pEnd)   *pEnd   = nC;
196   } else if (height == 1) {
197     /* Faces */
198     if (pStart) *pStart = nC+nV;
199     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
200   } else if (height == dim) {
201     /* Vertices */
202     if (pStart) *pStart = nC;
203     if (pEnd)   *pEnd   = nC+nV;
204   } else if (height < 0) {
205     /* All points */
206     if (pStart) *pStart = 0;
207     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
208   } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
209   PetscFunctionReturn(0);
210 }
211 
212 #undef __FUNCT__
213 #define __FUNCT__ "DMDAGetDepthStratum"
214 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
215 {
216   DM_DA         *da  = (DM_DA*) dm->data;
217   const PetscInt dim = da->dim;
218   PetscInt       nC, nV, nXF, nYF, nZF;
219   PetscErrorCode ierr;
220 
221   PetscFunctionBegin;
222   if (pStart) PetscValidIntPointer(pStart,3);
223   if (pEnd)   PetscValidIntPointer(pEnd,4);
224   ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
225   ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr);
226   ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr);
227   if (depth == dim) {
228     /* Cells */
229     if (pStart) *pStart = 0;
230     if (pEnd)   *pEnd   = nC;
231   } else if (depth == dim-1) {
232     /* Faces */
233     if (pStart) *pStart = nC+nV;
234     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
235   } else if (depth == 0) {
236     /* Vertices */
237     if (pStart) *pStart = nC;
238     if (pEnd)   *pEnd   = nC+nV;
239   } else if (depth < 0) {
240     /* All points */
241     if (pStart) *pStart = 0;
242     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
243   } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
244   PetscFunctionReturn(0);
245 }
246 
247 #undef __FUNCT__
248 #define __FUNCT__ "DMDAGetConeSize"
249 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
250 {
251   DM_DA         *da  = (DM_DA*) dm->data;
252   const PetscInt dim = da->dim;
253   PetscInt       nC, nV, nXF, nYF, nZF;
254   PetscErrorCode ierr;
255 
256   PetscFunctionBegin;
257   *coneSize = 0;
258   ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
259   ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr);
260   ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr);
261   switch (dim) {
262   case 2:
263     if (p >= 0) {
264       if (p < nC) {
265         *coneSize = 4;
266       } else if (p < nC+nV) {
267         *coneSize = 0;
268       } else if (p < nC+nV+nXF+nYF+nZF) {
269         *coneSize = 2;
270       } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
271     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
272     break;
273   case 3:
274     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
275     break;
276   }
277   PetscFunctionReturn(0);
278 }
279 
280 #undef __FUNCT__
281 #define __FUNCT__ "DMDAGetCone"
282 PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
283 {
284   DM_DA         *da  = (DM_DA*) dm->data;
285   const PetscInt dim = da->dim;
286   PetscInt       nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
287   PetscErrorCode ierr;
288 
289   PetscFunctionBegin;
290   if (!cone) {ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr);}
291   ierr = DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);CHKERRQ(ierr);
292   ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr);
293   ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr);
294   switch (dim) {
295   case 2:
296     if (p >= 0) {
297       if (p < nC) {
298         const PetscInt cy = p / nCx;
299         const PetscInt cx = p % nCx;
300 
301         (*cone)[0] = cy*nxF + cx + nC+nV;
302         (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
303         (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
304         (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
305         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
306       } else if (p < nC+nV) {
307       } else if (p < nC+nV+nXF) {
308         const PetscInt fy = (p - nC+nV) / nxF;
309         const PetscInt fx = (p - nC+nV) % nxF;
310 
311         (*cone)[0] = fy*nVx + fx + nC;
312         (*cone)[1] = fy*nVx + fx + 1 + nC;
313       } else if (p < nC+nV+nXF+nYF) {
314         const PetscInt fx = (p - nC+nV+nXF) / nyF;
315         const PetscInt fy = (p - nC+nV+nXF) % nyF;
316 
317         (*cone)[0] = fy*nVx + fx + nC;
318         (*cone)[1] = fy*nVx + fx + nVx + nC;
319       } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
320     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
321     break;
322   case 3:
323     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
324     break;
325   }
326   PetscFunctionReturn(0);
327 }
328 
329 #undef __FUNCT__
330 #define __FUNCT__ "DMDARestoreCone"
331 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
332 {
333   PetscErrorCode ierr;
334 
335   PetscFunctionBegin;
336   ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr);
337   PetscFunctionReturn(0);
338 }
339 
340 #undef __FUNCT__
341 #define __FUNCT__ "DMDACreateSection"
342 /*@C
343   DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with
344   different numbers of dofs on vertices, cells, and faces in each direction.
345 
346   Input Parameters:
347 + dm- The DMDA
348 . numFields - The number of fields
349 . numComp - The number of components in each field, or NULL for 1
350 . numVertexDof - The number of dofs per vertex for each field, or NULL
351 . numFaceDof - The number of dofs per face for each field and direction, or NULL
352 - numCellDof - The number of dofs per cell for each field, or NULL
353 
354   Level: developer
355 
356   Note:
357   The default DMDA numbering is as follows:
358 
359     - Cells:    [0,             nC)
360     - Vertices: [nC,            nC+nV)
361     - X-Faces:  [nC+nV,         nC+nV+nXF)         normal is +- x-dir
362     - Y-Faces:  [nC+nV+nXF,     nC+nV+nXF+nYF)     normal is +- y-dir
363     - Z-Faces:  [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir
364 
365   We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment.
366 @*/
367 PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numVertexDof[], PetscInt numFaceDof[], PetscInt numCellDof[])
368 {
369   DM_DA            *da  = (DM_DA*) dm->data;
370   PetscSection      section;
371   const PetscInt    dim = da->dim;
372   PetscInt          numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0};
373   PetscSF           sf;
374   PetscMPIInt       rank;
375   const PetscMPIInt *neighbors;
376   PetscInt          *localPoints;
377   PetscSFNode       *remotePoints;
378   PetscInt          nleaves = 0,  nleavesCheck = 0, nL = 0;
379   PetscInt          nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF;
380   PetscInt          pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd;
381   PetscInt          f, v, c, xf, yf, zf, xn, yn, zn;
382   PetscErrorCode    ierr;
383 
384   PetscFunctionBegin;
385   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
386   PetscValidPointer(s, 4);
387   ierr    = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
388   ierr    = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
389   ierr    = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr);
390   ierr    = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr);
391   ierr    = DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);CHKERRQ(ierr);
392   ierr    = DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);CHKERRQ(ierr);
393   ierr    = DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);CHKERRQ(ierr);
394   ierr    = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr);
395   xfStart = vEnd;  xfEnd = xfStart+nXF;
396   yfStart = xfEnd; yfEnd = yfStart+nYF;
397   zfStart = yfEnd; zfEnd = zfStart+nZF;
398   if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd);
399   /* Create local section */
400   ierr = DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);CHKERRQ(ierr);
401   for (f = 0; f < numFields; ++f) {
402     if (numVertexDof) numVertexTotDof  += numVertexDof[f];
403     if (numCellDof)   numCellTotDof    += numCellDof[f];
404     if (numFaceDof) {
405       numFaceTotDof[0] += numFaceDof[f*dim+0];
406       numFaceTotDof[1] += dim > 1 ? numFaceDof[f*dim+1] : 0;
407       numFaceTotDof[2] += dim > 2 ? numFaceDof[f*dim+2] : 0;
408     }
409   }
410   ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);CHKERRQ(ierr);
411   if (numFields > 1) {
412     ierr = PetscSectionSetNumFields(section, numFields);CHKERRQ(ierr);
413     for (f = 0; f < numFields; ++f) {
414       const char *name;
415 
416       ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr);
417       ierr = PetscSectionSetFieldName(section, f, name);CHKERRQ(ierr);
418       if (numComp) {
419         ierr = PetscSectionSetFieldComponents(section, f, numComp[f]);CHKERRQ(ierr);
420       }
421     }
422   } else {
423     numFields = 0;
424   }
425   ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr);
426   if (numVertexDof) {
427     for (v = vStart; v < vEnd; ++v) {
428       for (f = 0; f < numFields; ++f) {
429         ierr = PetscSectionSetFieldDof(section, v, f, numVertexDof[f]);CHKERRQ(ierr);
430       }
431       ierr = PetscSectionSetDof(section, v, numVertexTotDof);CHKERRQ(ierr);
432     }
433   }
434   if (numFaceDof) {
435     for (xf = xfStart; xf < xfEnd; ++xf) {
436       for (f = 0; f < numFields; ++f) {
437         ierr = PetscSectionSetFieldDof(section, xf, f, numFaceDof[f*dim+0]);CHKERRQ(ierr);
438       }
439       ierr = PetscSectionSetDof(section, xf, numFaceTotDof[0]);CHKERRQ(ierr);
440     }
441     for (yf = yfStart; yf < yfEnd; ++yf) {
442       for (f = 0; f < numFields; ++f) {
443         ierr = PetscSectionSetFieldDof(section, yf, f, numFaceDof[f*dim+1]);CHKERRQ(ierr);
444       }
445       ierr = PetscSectionSetDof(section, yf, numFaceTotDof[1]);CHKERRQ(ierr);
446     }
447     for (zf = zfStart; zf < zfEnd; ++zf) {
448       for (f = 0; f < numFields; ++f) {
449         ierr = PetscSectionSetFieldDof(section, zf, f, numFaceDof[f*dim+2]);CHKERRQ(ierr);
450       }
451       ierr = PetscSectionSetDof(section, zf, numFaceTotDof[2]);CHKERRQ(ierr);
452     }
453   }
454   if (numCellDof) {
455     for (c = cStart; c < cEnd; ++c) {
456       for (f = 0; f < numFields; ++f) {
457         ierr = PetscSectionSetFieldDof(section, c, f, numCellDof[f]);CHKERRQ(ierr);
458       }
459       ierr = PetscSectionSetDof(section, c, numCellTotDof);CHKERRQ(ierr);
460     }
461   }
462   ierr = PetscSectionSetUp(section);CHKERRQ(ierr);
463   /* Create mesh point SF */
464   ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr);
465   for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
466     for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
467       for (xn = 0; xn < 3; ++xn) {
468         const PetscInt xp       = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
469         const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
470 
471         if (neighbor >= 0 && neighbor < rank) {
472           nleaves += (!xp ? nVx : 1) * (!yp ? nVy : 1) * (!zp ? nVz : 1); /* vertices */
473           if (xp && !yp && !zp) {
474             nleaves += nxF; /* x faces */
475           } else if (yp && !zp && !xp) {
476             nleaves += nyF; /* y faces */
477           } else if (zp && !xp && !yp) {
478             nleaves += nzF; /* z faces */
479           }
480         }
481       }
482     }
483   }
484   ierr = PetscMalloc2(nleaves,&localPoints,nleaves,&remotePoints);CHKERRQ(ierr);
485   for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
486     for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
487       for (xn = 0; xn < 3; ++xn) {
488         const PetscInt xp       = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
489         const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
490         PetscInt       xv, yv, zv;
491 
492         if (neighbor >= 0 && neighbor < rank) {
493           if (xp < 0) { /* left */
494             if (yp < 0) { /* bottom */
495               if (zp < 0) { /* back */
496                 const PetscInt localVertex  = (      0*nVy +     0)*nVx +     0 + nC;
497                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
498                 nleavesCheck += 1; /* left bottom back vertex */
499 
500                 localPoints[nL]        = localVertex;
501                 remotePoints[nL].rank  = neighbor;
502                 remotePoints[nL].index = remoteVertex;
503                 ++nL;
504               } else if (zp > 0) { /* front */
505                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx +     0 + nC;
506                 const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
507                 nleavesCheck += 1; /* left bottom front vertex */
508 
509                 localPoints[nL]        = localVertex;
510                 remotePoints[nL].rank  = neighbor;
511                 remotePoints[nL].index = remoteVertex;
512                 ++nL;
513               } else {
514                 nleavesCheck += nVz; /* left bottom vertices */
515                 for (zv = 0; zv < nVz; ++zv, ++nL) {
516                   const PetscInt localVertex  = (zv*nVy +     0)*nVx +     0 + nC;
517                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
518 
519                   localPoints[nL]        = localVertex;
520                   remotePoints[nL].rank  = neighbor;
521                   remotePoints[nL].index = remoteVertex;
522                 }
523               }
524             } else if (yp > 0) { /* top */
525               if (zp < 0) { /* back */
526                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx +     0 + nC;
527                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
528                 nleavesCheck += 1; /* left top back vertex */
529 
530                 localPoints[nL]        = localVertex;
531                 remotePoints[nL].rank  = neighbor;
532                 remotePoints[nL].index = remoteVertex;
533                 ++nL;
534               } else if (zp > 0) { /* front */
535                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC;
536                 const PetscInt remoteVertex = (      0*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
537                 nleavesCheck += 1; /* left top front vertex */
538 
539                 localPoints[nL]        = localVertex;
540                 remotePoints[nL].rank  = neighbor;
541                 remotePoints[nL].index = remoteVertex;
542                 ++nL;
543               } else {
544                 nleavesCheck += nVz; /* left top vertices */
545                 for (zv = 0; zv < nVz; ++zv, ++nL) {
546                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx +     0 + nC;
547                   const PetscInt remoteVertex = (zv*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
548 
549                   localPoints[nL]        = localVertex;
550                   remotePoints[nL].rank  = neighbor;
551                   remotePoints[nL].index = remoteVertex;
552                 }
553               }
554             } else {
555               if (zp < 0) { /* back */
556                 nleavesCheck += nVy; /* left back vertices */
557                 for (yv = 0; yv < nVy; ++yv, ++nL) {
558                   const PetscInt localVertex  = (      0*nVy + yv)*nVx +     0 + nC;
559                   const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
560 
561                   localPoints[nL]        = localVertex;
562                   remotePoints[nL].rank  = neighbor;
563                   remotePoints[nL].index = remoteVertex;
564                 }
565               } else if (zp > 0) { /* front */
566                 nleavesCheck += nVy; /* left front vertices */
567                 for (yv = 0; yv < nVy; ++yv, ++nL) {
568                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx +     0 + nC;
569                   const PetscInt remoteVertex = (      0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
570 
571                   localPoints[nL]        = localVertex;
572                   remotePoints[nL].rank  = neighbor;
573                   remotePoints[nL].index = remoteVertex;
574                 }
575               } else {
576                 nleavesCheck += nVy*nVz; /* left vertices */
577                 for (zv = 0; zv < nVz; ++zv) {
578                   for (yv = 0; yv < nVy; ++yv, ++nL) {
579                     const PetscInt localVertex  = (zv*nVy + yv)*nVx +     0 + nC;
580                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
581 
582                     localPoints[nL]        = localVertex;
583                     remotePoints[nL].rank  = neighbor;
584                     remotePoints[nL].index = remoteVertex;
585                   }
586                 }
587                 nleavesCheck += nxF;     /* left faces */
588                 for (xf = 0; xf < nxF; ++xf, ++nL) {
589                   /* THIS IS WRONG */
590                   const PetscInt localFace  = 0 + nC+nV;
591                   const PetscInt remoteFace = 0 + nC+nV;
592 
593                   localPoints[nL]        = localFace;
594                   remotePoints[nL].rank  = neighbor;
595                   remotePoints[nL].index = remoteFace;
596                 }
597               }
598             }
599           } else if (xp > 0) { /* right */
600             if (yp < 0) { /* bottom */
601               if (zp < 0) { /* back */
602                 const PetscInt localVertex  = (      0*nVy +     0)*nVx + nVx-1 + nC;
603                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
604                 nleavesCheck += 1; /* right bottom back vertex */
605 
606                 localPoints[nL]        = localVertex;
607                 remotePoints[nL].rank  = neighbor;
608                 remotePoints[nL].index = remoteVertex;
609                 ++nL;
610               } else if (zp > 0) { /* front */
611                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC;
612                 const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
613                 nleavesCheck += 1; /* right bottom front vertex */
614 
615                 localPoints[nL]        = localVertex;
616                 remotePoints[nL].rank  = neighbor;
617                 remotePoints[nL].index = remoteVertex;
618                 ++nL;
619               } else {
620                 nleavesCheck += nVz; /* right bottom vertices */
621                 for (zv = 0; zv < nVz; ++zv, ++nL) {
622                   const PetscInt localVertex  = (zv*nVy +     0)*nVx + nVx-1 + nC;
623                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
624 
625                   localPoints[nL]        = localVertex;
626                   remotePoints[nL].rank  = neighbor;
627                   remotePoints[nL].index = remoteVertex;
628                 }
629               }
630             } else if (yp > 0) { /* top */
631               if (zp < 0) { /* back */
632                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + nVx-1 + nC;
633                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
634                 nleavesCheck += 1; /* right top back vertex */
635 
636                 localPoints[nL]        = localVertex;
637                 remotePoints[nL].rank  = neighbor;
638                 remotePoints[nL].index = remoteVertex;
639                 ++nL;
640               } else if (zp > 0) { /* front */
641                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC;
642                 const PetscInt remoteVertex = (      0*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
643                 nleavesCheck += 1; /* right top front vertex */
644 
645                 localPoints[nL]        = localVertex;
646                 remotePoints[nL].rank  = neighbor;
647                 remotePoints[nL].index = remoteVertex;
648                 ++nL;
649               } else {
650                 nleavesCheck += nVz; /* right top vertices */
651                 for (zv = 0; zv < nVz; ++zv, ++nL) {
652                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + nVx-1 + nC;
653                   const PetscInt remoteVertex = (zv*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
654 
655                   localPoints[nL]        = localVertex;
656                   remotePoints[nL].rank  = neighbor;
657                   remotePoints[nL].index = remoteVertex;
658                 }
659               }
660             } else {
661               if (zp < 0) { /* back */
662                 nleavesCheck += nVy; /* right back vertices */
663                 for (yv = 0; yv < nVy; ++yv, ++nL) {
664                   const PetscInt localVertex  = (      0*nVy + yv)*nVx + nVx-1 + nC;
665                   const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
666 
667                   localPoints[nL]        = localVertex;
668                   remotePoints[nL].rank  = neighbor;
669                   remotePoints[nL].index = remoteVertex;
670                 }
671               } else if (zp > 0) { /* front */
672                 nleavesCheck += nVy; /* right front vertices */
673                 for (yv = 0; yv < nVy; ++yv, ++nL) {
674                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC;
675                   const PetscInt remoteVertex = (      0*nVy + yv)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
676 
677                   localPoints[nL]        = localVertex;
678                   remotePoints[nL].rank  = neighbor;
679                   remotePoints[nL].index = remoteVertex;
680                 }
681               } else {
682                 nleavesCheck += nVy*nVz; /* right vertices */
683                 for (zv = 0; zv < nVz; ++zv) {
684                   for (yv = 0; yv < nVy; ++yv, ++nL) {
685                     const PetscInt localVertex  = (zv*nVy + yv)*nVx + nVx-1 + nC;
686                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0     + nC; /* TODO: Correct this for neighbor sizes */
687 
688                     localPoints[nL]        = localVertex;
689                     remotePoints[nL].rank  = neighbor;
690                     remotePoints[nL].index = remoteVertex;
691                   }
692                 }
693                 nleavesCheck += nxF;     /* right faces */
694                 for (xf = 0; xf < nxF; ++xf, ++nL) {
695                   /* THIS IS WRONG */
696                   const PetscInt localFace  = 0 + nC+nV;
697                   const PetscInt remoteFace = 0 + nC+nV;
698 
699                   localPoints[nL]        = localFace;
700                   remotePoints[nL].rank  = neighbor;
701                   remotePoints[nL].index = remoteFace;
702                 }
703               }
704             }
705           } else {
706             if (yp < 0) { /* bottom */
707               if (zp < 0) { /* back */
708                 nleavesCheck += nVx; /* bottom back vertices */
709                 for (xv = 0; xv < nVx; ++xv, ++nL) {
710                   const PetscInt localVertex  = (      0*nVy +     0)*nVx + xv + nC;
711                   const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
712 
713                   localPoints[nL]        = localVertex;
714                   remotePoints[nL].rank  = neighbor;
715                   remotePoints[nL].index = remoteVertex;
716                 }
717               } else if (zp > 0) { /* front */
718                 nleavesCheck += nVx; /* bottom front vertices */
719                 for (xv = 0; xv < nVx; ++xv, ++nL) {
720                   const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + xv + nC;
721                   const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
722 
723                   localPoints[nL]        = localVertex;
724                   remotePoints[nL].rank  = neighbor;
725                   remotePoints[nL].index = remoteVertex;
726                 }
727               } else {
728                 nleavesCheck += nVx*nVz; /* bottom vertices */
729                 for (zv = 0; zv < nVz; ++zv) {
730                   for (xv = 0; xv < nVx; ++xv, ++nL) {
731                     const PetscInt localVertex  = (zv*nVy +     0)*nVx + xv + nC;
732                     const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
733 
734                     localPoints[nL]        = localVertex;
735                     remotePoints[nL].rank  = neighbor;
736                     remotePoints[nL].index = remoteVertex;
737                   }
738                 }
739                 nleavesCheck += nyF;     /* bottom faces */
740                 for (yf = 0; yf < nyF; ++yf, ++nL) {
741                   /* THIS IS WRONG */
742                   const PetscInt localFace  = 0 + nC+nV;
743                   const PetscInt remoteFace = 0 + nC+nV;
744 
745                   localPoints[nL]        = localFace;
746                   remotePoints[nL].rank  = neighbor;
747                   remotePoints[nL].index = remoteFace;
748                 }
749               }
750             } else if (yp > 0) { /* top */
751               if (zp < 0) { /* back */
752                 nleavesCheck += nVx; /* top back vertices */
753                 for (xv = 0; xv < nVx; ++xv, ++nL) {
754                   const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + xv + nC;
755                   const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
756 
757                   localPoints[nL]        = localVertex;
758                   remotePoints[nL].rank  = neighbor;
759                   remotePoints[nL].index = remoteVertex;
760                 }
761               } else if (zp > 0) { /* front */
762                 nleavesCheck += nVx; /* top front vertices */
763                 for (xv = 0; xv < nVx; ++xv, ++nL) {
764                   const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC;
765                   const PetscInt remoteVertex = (      0*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
766 
767                   localPoints[nL]        = localVertex;
768                   remotePoints[nL].rank  = neighbor;
769                   remotePoints[nL].index = remoteVertex;
770                 }
771               } else {
772                 nleavesCheck += nVx*nVz; /* top vertices */
773                 for (zv = 0; zv < nVz; ++zv) {
774                   for (xv = 0; xv < nVx; ++xv, ++nL) {
775                     const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + xv + nC;
776                     const PetscInt remoteVertex = (zv*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
777 
778                     localPoints[nL]        = localVertex;
779                     remotePoints[nL].rank  = neighbor;
780                     remotePoints[nL].index = remoteVertex;
781                   }
782                 }
783                 nleavesCheck += nyF;     /* top faces */
784                 for (yf = 0; yf < nyF; ++yf, ++nL) {
785                   /* THIS IS WRONG */
786                   const PetscInt localFace  = 0 + nC+nV;
787                   const PetscInt remoteFace = 0 + nC+nV;
788 
789                   localPoints[nL]        = localFace;
790                   remotePoints[nL].rank  = neighbor;
791                   remotePoints[nL].index = remoteFace;
792                 }
793               }
794             } else {
795               if (zp < 0) { /* back */
796                 nleavesCheck += nVx*nVy; /* back vertices */
797                 for (yv = 0; yv < nVy; ++yv) {
798                   for (xv = 0; xv < nVx; ++xv, ++nL) {
799                     const PetscInt localVertex  = (      0*nVy + yv)*nVx + xv + nC;
800                     const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
801 
802                     localPoints[nL]        = localVertex;
803                     remotePoints[nL].rank  = neighbor;
804                     remotePoints[nL].index = remoteVertex;
805                   }
806                 }
807                 nleavesCheck += nzF;     /* back faces */
808                 for (zf = 0; zf < nzF; ++zf, ++nL) {
809                   /* THIS IS WRONG */
810                   const PetscInt localFace  = 0 + nC+nV;
811                   const PetscInt remoteFace = 0 + nC+nV;
812 
813                   localPoints[nL]        = localFace;
814                   remotePoints[nL].rank  = neighbor;
815                   remotePoints[nL].index = remoteFace;
816                 }
817               } else if (zp > 0) { /* front */
818                 nleavesCheck += nVx*nVy; /* front vertices */
819                 for (yv = 0; yv < nVy; ++yv) {
820                   for (xv = 0; xv < nVx; ++xv, ++nL) {
821                     const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + xv + nC;
822                     const PetscInt remoteVertex = (      0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
823 
824                     localPoints[nL]        = localVertex;
825                     remotePoints[nL].rank  = neighbor;
826                     remotePoints[nL].index = remoteVertex;
827                   }
828                 }
829                 nleavesCheck += nzF;     /* front faces */
830                 for (zf = 0; zf < nzF; ++zf, ++nL) {
831                   /* THIS IS WRONG */
832                   const PetscInt localFace  = 0 + nC+nV;
833                   const PetscInt remoteFace = 0 + nC+nV;
834 
835                   localPoints[nL]        = localFace;
836                   remotePoints[nL].rank  = neighbor;
837                   remotePoints[nL].index = remoteFace;
838                 }
839               } else {
840                 /* Nothing is shared from the interior */
841               }
842             }
843           }
844         }
845       }
846     }
847   }
848   /* TODO: Remove duplication in leaf determination */
849   if (nleaves != nleavesCheck) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck);
850   ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);CHKERRQ(ierr);
851   ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr);
852   ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr);
853   ierr = PetscSFDestroy(&sf);CHKERRQ(ierr);
854   ierr = DMSetDefaultSection(dm, section);CHKERRQ(ierr);
855 #undef __FUNCT__
856 #define __FUNCT__ "DMDASetVertexCoordinates"
857 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
858 {
859   DM_DA         *da = (DM_DA *) dm->data;
860   Vec            coordinates;
861   PetscSection   section;
862   PetscScalar   *coords;
863   PetscReal      h[3];
864   PetscInt       dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
865   PetscErrorCode ierr;
866 
867   PetscFunctionBegin;
868   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
869   ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
870   h[0] = (xu - xl)/M;
871   h[1] = (yu - yl)/N;
872   h[2] = (zu - zl)/P;
873   ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
874   ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr);
875   ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);CHKERRQ(ierr);
876   ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr);
877   ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr);
878   ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr);
879   for (v = vStart; v < vEnd; ++v) {
880     ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr);
881   }
882   ierr = PetscSectionSetUp(section);CHKERRQ(ierr);
883   ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
884   ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr);
885   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
886   for (k = 0; k < nVz; ++k) {
887     PetscInt ind[3] = {0, 0, k + da->zs}, d, off;
888 
889     for (j = 0; j < nVy; ++j) {
890       ind[1] = j + da->ys;
891       for (i = 0; i < nVx; ++i) {
892         const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;
893 
894         ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr);
895         ind[0] = i + da->xs;
896         for (d = 0; d < dim; ++d) {
897           coords[off+d] = h[d]*ind[d];
898         }
899       }
900     }
901   }
902   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
903   ierr = DMSetCoordinateSection(dm, section);CHKERRQ(ierr);
904   ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
905   ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);
906   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
907   PetscFunctionReturn(0);
908 }
909   PetscFunctionReturn(0);
910 }
911 
912 /* ------------------------------------------------------------------- */
913 
914 #undef __FUNCT__
915 #define __FUNCT__ "DMDAGetArray"
916 /*@C
917      DMDAGetArray - Gets a work array for a DMDA
918 
919     Input Parameter:
920 +    da - information about my local patch
921 -    ghosted - do you want arrays for the ghosted or nonghosted patch
922 
923     Output Parameters:
924 .    vptr - array data structured
925 
926     Note:  The vector values are NOT initialized and may have garbage in them, so you may need
927            to zero them.
928 
929   Level: advanced
930 
931 .seealso: DMDARestoreArray()
932 
933 @*/
934 PetscErrorCode  DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
935 {
936   PetscErrorCode ierr;
937   PetscInt       j,i,xs,ys,xm,ym,zs,zm;
938   char           *iarray_start;
939   void           **iptr = (void**)vptr;
940   DM_DA          *dd    = (DM_DA*)da->data;
941 
942   PetscFunctionBegin;
943   PetscValidHeaderSpecific(da,DM_CLASSID,1);
944   if (ghosted) {
945     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
946       if (dd->arrayghostedin[i]) {
947         *iptr                 = dd->arrayghostedin[i];
948         iarray_start          = (char*)dd->startghostedin[i];
949         dd->arrayghostedin[i] = NULL;
950         dd->startghostedin[i] = NULL;
951 
952         goto done;
953       }
954     }
955     xs = dd->Xs;
956     ys = dd->Ys;
957     zs = dd->Zs;
958     xm = dd->Xe-dd->Xs;
959     ym = dd->Ye-dd->Ys;
960     zm = dd->Ze-dd->Zs;
961   } else {
962     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
963       if (dd->arrayin[i]) {
964         *iptr          = dd->arrayin[i];
965         iarray_start   = (char*)dd->startin[i];
966         dd->arrayin[i] = NULL;
967         dd->startin[i] = NULL;
968 
969         goto done;
970       }
971     }
972     xs = dd->xs;
973     ys = dd->ys;
974     zs = dd->zs;
975     xm = dd->xe-dd->xs;
976     ym = dd->ye-dd->ys;
977     zm = dd->ze-dd->zs;
978   }
979 
980   switch (dd->dim) {
981   case 1: {
982     void *ptr;
983 
984     ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
985 
986     ptr   = (void*)(iarray_start - xs*sizeof(PetscScalar));
987     *iptr = (void*)ptr;
988     break;
989   }
990   case 2: {
991     void **ptr;
992 
993     ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
994 
995     ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
996     for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
997     *iptr = (void*)ptr;
998     break;
999   }
1000   case 3: {
1001     void ***ptr,**bptr;
1002 
1003     ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
1004 
1005     ptr  = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1006     bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1007     for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1008     for (i=zs; i<zs+zm; i++) {
1009       for (j=ys; j<ys+ym; j++) {
1010         ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1011       }
1012     }
1013     *iptr = (void*)ptr;
1014     break;
1015   }
1016   default:
1017     SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
1018   }
1019 
1020 done:
1021   /* add arrays to the checked out list */
1022   if (ghosted) {
1023     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1024       if (!dd->arrayghostedout[i]) {
1025         dd->arrayghostedout[i] = *iptr;
1026         dd->startghostedout[i] = iarray_start;
1027         break;
1028       }
1029     }
1030   } else {
1031     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1032       if (!dd->arrayout[i]) {
1033         dd->arrayout[i] = *iptr;
1034         dd->startout[i] = iarray_start;
1035         break;
1036       }
1037     }
1038   }
1039   PetscFunctionReturn(0);
1040 }
1041 
1042 #undef __FUNCT__
1043 #define __FUNCT__ "DMDARestoreArray"
1044 /*@C
1045      DMDARestoreArray - Restores an array of derivative types for a DMDA
1046 
1047     Input Parameter:
1048 +    da - information about my local patch
1049 .    ghosted - do you want arrays for the ghosted or nonghosted patch
1050 -    vptr - array data structured to be passed to ad_FormFunctionLocal()
1051 
1052      Level: advanced
1053 
1054 .seealso: DMDAGetArray()
1055 
1056 @*/
1057 PetscErrorCode  DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1058 {
1059   PetscInt i;
1060   void     **iptr = (void**)vptr,*iarray_start = 0;
1061   DM_DA    *dd    = (DM_DA*)da->data;
1062 
1063   PetscFunctionBegin;
1064   PetscValidHeaderSpecific(da,DM_CLASSID,1);
1065   if (ghosted) {
1066     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1067       if (dd->arrayghostedout[i] == *iptr) {
1068         iarray_start           = dd->startghostedout[i];
1069         dd->arrayghostedout[i] = NULL;
1070         dd->startghostedout[i] = NULL;
1071         break;
1072       }
1073     }
1074     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1075       if (!dd->arrayghostedin[i]) {
1076         dd->arrayghostedin[i] = *iptr;
1077         dd->startghostedin[i] = iarray_start;
1078         break;
1079       }
1080     }
1081   } else {
1082     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1083       if (dd->arrayout[i] == *iptr) {
1084         iarray_start    = dd->startout[i];
1085         dd->arrayout[i] = NULL;
1086         dd->startout[i] = NULL;
1087         break;
1088       }
1089     }
1090     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1091       if (!dd->arrayin[i]) {
1092         dd->arrayin[i] = *iptr;
1093         dd->startin[i] = iarray_start;
1094         break;
1095       }
1096     }
1097   }
1098   PetscFunctionReturn(0);
1099 }
1100 
1101