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