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