xref: /petsc/src/dm/impls/da/dalocal.c (revision d089267096544a113110f12e19de0de1dfbc8169) 
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__ "DMDAGetDepthStratum"
215 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
216 {
217   DM_DA         *da  = (DM_DA*) dm->data;
218   const PetscInt dim = da->dim;
219   PetscInt       nC, nV, nXF, nYF, nZF;
220   PetscErrorCode ierr;
221 
222   PetscFunctionBegin;
223   if (pStart) PetscValidIntPointer(pStart,3);
224   if (pEnd)   PetscValidIntPointer(pEnd,4);
225   ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
226   ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr);
227   ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr);
228   if (depth == dim) {
229     /* Cells */
230     if (pStart) *pStart = 0;
231     if (pEnd)   *pEnd   = nC;
232   } else if (depth == dim-1) {
233     /* Faces */
234     if (pStart) *pStart = nC+nV;
235     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
236   } else if (depth == 0) {
237     /* Vertices */
238     if (pStart) *pStart = nC;
239     if (pEnd)   *pEnd   = nC+nV;
240   } else if (depth < 0) {
241     /* All points */
242     if (pStart) *pStart = 0;
243     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
244   } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
245   PetscFunctionReturn(0);
246 }
247 
248 #undef __FUNCT__
249 #define __FUNCT__ "DMDAGetConeSize"
250 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
251 {
252   DM_DA         *da  = (DM_DA*) dm->data;
253   const PetscInt dim = da->dim;
254   PetscInt       nC, nV, nXF, nYF, nZF;
255   PetscErrorCode ierr;
256 
257   PetscFunctionBegin;
258   *coneSize = 0;
259   ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
260   ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr);
261   ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr);
262   switch (dim) {
263   case 2:
264     if (p >= 0) {
265       if (p < nC) {
266         *coneSize = 4;
267       } else if (p < nC+nV) {
268         *coneSize = 0;
269       } else if (p < nC+nV+nXF+nYF+nZF) {
270         *coneSize = 2;
271       } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
272     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
273     break;
274   case 3:
275     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
276     break;
277   }
278   PetscFunctionReturn(0);
279 }
280 
281 #undef __FUNCT__
282 #define __FUNCT__ "DMDAGetCone"
283 PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
284 {
285   DM_DA         *da  = (DM_DA*) dm->data;
286   const PetscInt dim = da->dim;
287   PetscInt       nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
288   PetscErrorCode ierr;
289 
290   PetscFunctionBegin;
291   if (!cone) {ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr);}
292   ierr = DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);CHKERRQ(ierr);
293   ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr);
294   ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr);
295   switch (dim) {
296   case 2:
297     if (p >= 0) {
298       if (p < nC) {
299         const PetscInt cy = p / nCx;
300         const PetscInt cx = p % nCx;
301 
302         (*cone)[0] = cy*nxF + cx + nC+nV;
303         (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
304         (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
305         (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
306         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
307       } else if (p < nC+nV) {
308       } else if (p < nC+nV+nXF) {
309         const PetscInt fy = (p - nC+nV) / nxF;
310         const PetscInt fx = (p - nC+nV) % nxF;
311 
312         (*cone)[0] = fy*nVx + fx + nC;
313         (*cone)[1] = fy*nVx + fx + 1 + nC;
314       } else if (p < nC+nV+nXF+nYF) {
315         const PetscInt fx = (p - nC+nV+nXF) / nyF;
316         const PetscInt fy = (p - nC+nV+nXF) % nyF;
317 
318         (*cone)[0] = fy*nVx + fx + nC;
319         (*cone)[1] = fy*nVx + fx + nVx + nC;
320       } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
321     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
322     break;
323   case 3:
324     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
325     break;
326   }
327   PetscFunctionReturn(0);
328 }
329 
330 #undef __FUNCT__
331 #define __FUNCT__ "DMDARestoreCone"
332 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
333 {
334   PetscErrorCode ierr;
335 
336   PetscFunctionBegin;
337   ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr);
338   PetscFunctionReturn(0);
339 }
340 
341 #undef __FUNCT__
342 #define __FUNCT__ "DMDACreateSection"
343 /*@C
344   DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with
345   different numbers of dofs on vertices, cells, and faces in each direction.
346 
347   Input Parameters:
348 + dm- The DMDA
349 . numFields - The number of fields
350 . numComp - The number of components in each field, or NULL for 1
351 . numVertexDof - The number of dofs per vertex for each field, or NULL
352 . numFaceDof - The number of dofs per face for each field and direction, or NULL
353 - numCellDof - The number of dofs per cell for each field, or NULL
354 
355   Level: developer
356 
357   Note:
358   The default DMDA numbering is as follows:
359 
360     - Cells:    [0,             nC)
361     - Vertices: [nC,            nC+nV)
362     - X-Faces:  [nC+nV,         nC+nV+nXF)         normal is +- x-dir
363     - Y-Faces:  [nC+nV+nXF,     nC+nV+nXF+nYF)     normal is +- y-dir
364     - Z-Faces:  [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir
365 
366   We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment.
367 @*/
368 PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numVertexDof[], PetscInt numFaceDof[], PetscInt numCellDof[])
369 {
370   DM_DA            *da  = (DM_DA*) dm->data;
371   PetscSection      section;
372   const PetscInt    dim = da->dim;
373   PetscInt          numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0};
374   PetscBT           isLeaf;
375   PetscSF           sf;
376   PetscMPIInt       rank;
377   const PetscMPIInt *neighbors;
378   PetscInt          *localPoints;
379   PetscSFNode       *remotePoints;
380   PetscInt          nleaves = 0,  nleavesCheck = 0, nL = 0;
381   PetscInt          nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF;
382   PetscInt          pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd;
383   PetscInt          f, v, c, xf, yf, zf, xn, yn, zn;
384   PetscErrorCode    ierr;
385 
386   PetscFunctionBegin;
387   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
388   PetscValidPointer(s, 4);
389   ierr    = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
390   ierr    = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr);
391   ierr    = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr);
392   ierr    = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr);
393   ierr    = DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);CHKERRQ(ierr);
394   ierr    = DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);CHKERRQ(ierr);
395   ierr    = DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);CHKERRQ(ierr);
396   ierr    = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr);
397   xfStart = vEnd;  xfEnd = xfStart+nXF;
398   yfStart = xfEnd; yfEnd = yfStart+nYF;
399   zfStart = yfEnd; zfEnd = zfStart+nZF;
400   if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd);
401   /* Create local section */
402   ierr = DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);CHKERRQ(ierr);
403   for (f = 0; f < numFields; ++f) {
404     if (numVertexDof) numVertexTotDof  += numVertexDof[f];
405     if (numCellDof)   numCellTotDof    += numCellDof[f];
406     if (numFaceDof) {
407       numFaceTotDof[0] += numFaceDof[f*dim+0];
408       numFaceTotDof[1] += dim > 1 ? numFaceDof[f*dim+1] : 0;
409       numFaceTotDof[2] += dim > 2 ? numFaceDof[f*dim+2] : 0;
410     }
411   }
412   ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);CHKERRQ(ierr);
413   if (numFields > 1) {
414     ierr = PetscSectionSetNumFields(section, numFields);CHKERRQ(ierr);
415     for (f = 0; f < numFields; ++f) {
416       const char *name;
417 
418       ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr);
419       ierr = PetscSectionSetFieldName(section, f, name);CHKERRQ(ierr);
420       if (numComp) {
421         ierr = PetscSectionSetFieldComponents(section, f, numComp[f]);CHKERRQ(ierr);
422       }
423     }
424   } else {
425     numFields = 0;
426   }
427   ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr);
428   if (numVertexDof) {
429     for (v = vStart; v < vEnd; ++v) {
430       for (f = 0; f < numFields; ++f) {
431         ierr = PetscSectionSetFieldDof(section, v, f, numVertexDof[f]);CHKERRQ(ierr);
432       }
433       ierr = PetscSectionSetDof(section, v, numVertexTotDof);CHKERRQ(ierr);
434     }
435   }
436   if (numFaceDof) {
437     for (xf = xfStart; xf < xfEnd; ++xf) {
438       for (f = 0; f < numFields; ++f) {
439         ierr = PetscSectionSetFieldDof(section, xf, f, numFaceDof[f*dim+0]);CHKERRQ(ierr);
440       }
441       ierr = PetscSectionSetDof(section, xf, numFaceTotDof[0]);CHKERRQ(ierr);
442     }
443     for (yf = yfStart; yf < yfEnd; ++yf) {
444       for (f = 0; f < numFields; ++f) {
445         ierr = PetscSectionSetFieldDof(section, yf, f, numFaceDof[f*dim+1]);CHKERRQ(ierr);
446       }
447       ierr = PetscSectionSetDof(section, yf, numFaceTotDof[1]);CHKERRQ(ierr);
448     }
449     for (zf = zfStart; zf < zfEnd; ++zf) {
450       for (f = 0; f < numFields; ++f) {
451         ierr = PetscSectionSetFieldDof(section, zf, f, numFaceDof[f*dim+2]);CHKERRQ(ierr);
452       }
453       ierr = PetscSectionSetDof(section, zf, numFaceTotDof[2]);CHKERRQ(ierr);
454     }
455   }
456   if (numCellDof) {
457     for (c = cStart; c < cEnd; ++c) {
458       for (f = 0; f < numFields; ++f) {
459         ierr = PetscSectionSetFieldDof(section, c, f, numCellDof[f]);CHKERRQ(ierr);
460       }
461       ierr = PetscSectionSetDof(section, c, numCellTotDof);CHKERRQ(ierr);
462     }
463   }
464   ierr = PetscSectionSetUp(section);CHKERRQ(ierr);
465   /* Create mesh point SF */
466   ierr = PetscBTCreate(pEnd-pStart, &isLeaf);CHKERRQ(ierr);
467   ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr);
468   for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
469     for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
470       for (xn = 0; xn < 3; ++xn) {
471         const PetscInt xp       = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
472         const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
473         PetscInt       xv, yv, zv;
474 
475         if (neighbor >= 0 && neighbor < rank) {
476           if (xp < 0) { /* left */
477             if (yp < 0) { /* bottom */
478               if (zp < 0) { /* back */
479                 const PetscInt localVertex  = (      0*nVy +     0)*nVx +     0 + nC; /* left bottom back vertex */
480                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
481               } else if (zp > 0) { /* front */
482                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* left bottom front vertex */
483                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
484               } else {
485                 for (zv = 0; zv < nVz; ++zv) {
486                   const PetscInt localVertex  = (zv*nVy +     0)*nVx +     0 + nC; /* left bottom vertices */
487                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
488                 }
489               }
490             } else if (yp > 0) { /* top */
491               if (zp < 0) { /* back */
492                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx +     0 + nC; /* left top back vertex */
493                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
494               } else if (zp > 0) { /* front */
495                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* left top front vertex */
496                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
497               } else {
498                 for (zv = 0; zv < nVz; ++zv) {
499                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx +     0 + nC; /* left top vertices */
500                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
501                 }
502               }
503             } else {
504               if (zp < 0) { /* back */
505                 for (yv = 0; yv < nVy; ++yv) {
506                   const PetscInt localVertex  = (      0*nVy + yv)*nVx +     0 + nC; /* left back vertices */
507                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
508                 }
509               } else if (zp > 0) { /* front */
510                 for (yv = 0; yv < nVy; ++yv) {
511                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx +     0 + nC; /* left front vertices */
512                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
513                 }
514               } else {
515                 for (zv = 0; zv < nVz; ++zv) {
516                   for (yv = 0; yv < nVy; ++yv) {
517                     const PetscInt localVertex  = (zv*nVy + yv)*nVx +     0 + nC; /* left vertices */
518                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
519                   }
520                 }
521 #if 0
522                 for (xf = 0; xf < nxF; ++xf) {
523                   /* THIS IS WRONG */
524                   const PetscInt localFace  = 0 + nC+nV; /* left faces */
525                   if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
526                 }
527 #endif
528               }
529             }
530           } else if (xp > 0) { /* right */
531             if (yp < 0) { /* bottom */
532               if (zp < 0) { /* back */
533                 const PetscInt localVertex  = (      0*nVy +     0)*nVx + nVx-1 + nC; /* right bottom back vertex */
534                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
535               } else if (zp > 0) { /* front */
536                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* right bottom front vertex */
537                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
538               } else {
539                 for (zv = 0; zv < nVz; ++zv) {
540                   const PetscInt localVertex  = (zv*nVy +     0)*nVx + nVx-1 + nC; /* right bottom vertices */
541                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
542                 }
543               }
544             } else if (yp > 0) { /* top */
545               if (zp < 0) { /* back */
546                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
547                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
548               } else if (zp > 0) { /* front */
549                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
550                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
551               } else {
552                 for (zv = 0; zv < nVz; ++zv) {
553                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
554                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
555                 }
556               }
557             } else {
558               if (zp < 0) { /* back */
559                 for (yv = 0; yv < nVy; ++yv) {
560                   const PetscInt localVertex  = (      0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
561                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
562                 }
563               } else if (zp > 0) { /* front */
564                 for (yv = 0; yv < nVy; ++yv) {
565                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
566                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
567                 }
568               } else {
569                 for (zv = 0; zv < nVz; ++zv) {
570                   for (yv = 0; yv < nVy; ++yv) {
571                     const PetscInt localVertex  = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
572                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
573                   }
574                 }
575 #if 0
576                 for (xf = 0; xf < nxF; ++xf) {
577                   /* THIS IS WRONG */
578                   const PetscInt localFace  = 0 + nC+nV; /* right faces */
579                   if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
580                 }
581 #endif
582               }
583             }
584           } else {
585             if (yp < 0) { /* bottom */
586               if (zp < 0) { /* back */
587                 for (xv = 0; xv < nVx; ++xv) {
588                   const PetscInt localVertex  = (      0*nVy +     0)*nVx + xv + nC; /* bottom back vertices */
589                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
590                 }
591               } else if (zp > 0) { /* front */
592                 for (xv = 0; xv < nVx; ++xv) {
593                   const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* bottom front vertices */
594                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
595                 }
596               } else {
597                 for (zv = 0; zv < nVz; ++zv) {
598                   for (xv = 0; xv < nVx; ++xv) {
599                     const PetscInt localVertex  = (zv*nVy +     0)*nVx + xv + nC; /* bottom vertices */
600                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
601                   }
602                 }
603 #if 0
604                 for (yf = 0; yf < nyF; ++yf) {
605                   /* THIS IS WRONG */
606                   const PetscInt localFace  = 0 + nC+nV; /* bottom faces */
607                   if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves;
608                 }
609 #endif
610               }
611             } else if (yp > 0) { /* top */
612               if (zp < 0) { /* back */
613                 for (xv = 0; xv < nVx; ++xv) {
614                   const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
615                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
616                 }
617               } else if (zp > 0) { /* front */
618                 for (xv = 0; xv < nVx; ++xv) {
619                   const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
620                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
621                 }
622               } else {
623                 for (zv = 0; zv < nVz; ++zv) {
624                   for (xv = 0; xv < nVx; ++xv) {
625                     const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
626                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
627                   }
628                 }
629 #if 0
630                 for (yf = 0; yf < nyF; ++yf) {
631                   /* THIS IS WRONG */
632                   const PetscInt localFace  = 0 + nC+nV; /* top faces */
633                   if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves;
634                 }
635 #endif
636               }
637             } else {
638               if (zp < 0) { /* back */
639                 for (yv = 0; yv < nVy; ++yv) {
640                   for (xv = 0; xv < nVx; ++xv) {
641                     const PetscInt localVertex  = (      0*nVy + yv)*nVx + xv + nC; /* back vertices */
642                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
643                   }
644                 }
645 #if 0
646                 for (zf = 0; zf < nzF; ++zf) {
647                   /* THIS IS WRONG */
648                   const PetscInt localFace  = 0 + nC+nV; /* back faces */
649                   if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
650                 }
651 #endif
652               } else if (zp > 0) { /* front */
653                 for (yv = 0; yv < nVy; ++yv) {
654                   for (xv = 0; xv < nVx; ++xv) {
655                     const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
656                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
657                   }
658                 }
659 #if 0
660                 for (zf = 0; zf < nzF; ++zf) {
661                   /* THIS IS WRONG */
662                   const PetscInt localFace  = 0 + nC+nV; /* front faces */
663                   if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
664                 }
665 #endif
666               } else {
667                 /* Nothing is shared from the interior */
668               }
669             }
670           }
671         }
672       }
673     }
674   }
675   ierr = PetscBTMemzero(pEnd-pStart, isLeaf);CHKERRQ(ierr);
676   ierr = PetscMalloc2(nleaves,PetscInt,&localPoints,nleaves,PetscSFNode,&remotePoints);CHKERRQ(ierr);
677   for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
678     for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
679       for (xn = 0; xn < 3; ++xn) {
680         const PetscInt xp       = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
681         const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
682         PetscInt       xv, yv, zv;
683 
684         if (neighbor >= 0 && neighbor < rank) {
685           if (xp < 0) { /* left */
686             if (yp < 0) { /* bottom */
687               if (zp < 0) { /* back */
688                 const PetscInt localVertex  = (      0*nVy +     0)*nVx +     0 + nC; /* left bottom back vertex */
689                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
690 
691                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
692                   localPoints[nL]        = localVertex;
693                   remotePoints[nL].rank  = neighbor;
694                   remotePoints[nL].index = remoteVertex;
695                   ++nL;
696                 }
697               } else if (zp > 0) { /* front */
698                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* left bottom front vertex */
699                 const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
700 
701                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
702                   localPoints[nL]        = localVertex;
703                   remotePoints[nL].rank  = neighbor;
704                   remotePoints[nL].index = remoteVertex;
705                   ++nL;
706                 }
707               } else {
708                 for (zv = 0; zv < nVz; ++zv) {
709                   const PetscInt localVertex  = (zv*nVy +     0)*nVx +     0 + nC; /* left bottom vertices */
710                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
711 
712                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
713                     localPoints[nL]        = localVertex;
714                     remotePoints[nL].rank  = neighbor;
715                     remotePoints[nL].index = remoteVertex;
716                     ++nL;
717                   }
718                 }
719               }
720             } else if (yp > 0) { /* top */
721               if (zp < 0) { /* back */
722                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx +     0 + nC; /* left top back vertex */
723                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
724 
725                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
726                   localPoints[nL]        = localVertex;
727                   remotePoints[nL].rank  = neighbor;
728                   remotePoints[nL].index = remoteVertex;
729                   ++nL;
730                 }
731               } else if (zp > 0) { /* front */
732                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* left top front vertex */
733                 const PetscInt remoteVertex = (      0*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
734 
735                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
736                   localPoints[nL]        = localVertex;
737                   remotePoints[nL].rank  = neighbor;
738                   remotePoints[nL].index = remoteVertex;
739                   ++nL;
740                 }
741               } else {
742                 for (zv = 0; zv < nVz; ++zv) {
743                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx +     0 + nC; /* left top vertices */
744                   const PetscInt remoteVertex = (zv*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
745 
746                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
747                     localPoints[nL]        = localVertex;
748                     remotePoints[nL].rank  = neighbor;
749                     remotePoints[nL].index = remoteVertex;
750                     ++nL;
751                   }
752                 }
753               }
754             } else {
755               if (zp < 0) { /* back */
756                 for (yv = 0; yv < nVy; ++yv) {
757                   const PetscInt localVertex  = (      0*nVy + yv)*nVx +     0 + nC; /* left back vertices */
758                   const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
759 
760                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
761                     localPoints[nL]        = localVertex;
762                     remotePoints[nL].rank  = neighbor;
763                     remotePoints[nL].index = remoteVertex;
764                     ++nL;
765                   }
766                 }
767               } else if (zp > 0) { /* front */
768                 for (yv = 0; yv < nVy; ++yv) {
769                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx +     0 + nC; /* left front vertices */
770                   const PetscInt remoteVertex = (      0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
771 
772                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
773                     localPoints[nL]        = localVertex;
774                     remotePoints[nL].rank  = neighbor;
775                     remotePoints[nL].index = remoteVertex;
776                     ++nL;
777                   }
778                 }
779               } else {
780                 for (zv = 0; zv < nVz; ++zv) {
781                   for (yv = 0; yv < nVy; ++yv) {
782                     const PetscInt localVertex  = (zv*nVy + yv)*nVx +     0 + nC; /* left vertices */
783                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
784 
785                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
786                       localPoints[nL]        = localVertex;
787                       remotePoints[nL].rank  = neighbor;
788                       remotePoints[nL].index = remoteVertex;
789                       ++nL;
790                     }
791                   }
792                 }
793 #if 0
794                 for (xf = 0; xf < nxF; ++xf) {
795                   /* THIS IS WRONG */
796                   const PetscInt localFace  = 0 + nC+nV; /* left faces */
797                   const PetscInt remoteFace = 0 + nC+nV;
798 
799                   if (!PetscBTLookupSet(isLeaf, localFace)) {
800                     localPoints[nL]        = localFace;
801                     remotePoints[nL].rank  = neighbor;
802                     remotePoints[nL].index = remoteFace;
803                   }
804                 }
805 #endif
806               }
807             }
808           } else if (xp > 0) { /* right */
809             if (yp < 0) { /* bottom */
810               if (zp < 0) { /* back */
811                 const PetscInt localVertex  = (      0*nVy +     0)*nVx + nVx-1 + nC; /* right bottom back vertex */
812                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
813 
814                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
815                   localPoints[nL]        = localVertex;
816                   remotePoints[nL].rank  = neighbor;
817                   remotePoints[nL].index = remoteVertex;
818                   ++nL;
819                 }
820               } else if (zp > 0) { /* front */
821                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* right bottom front vertex */
822                 const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx +     0 + 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               } else {
831                 nleavesCheck += nVz;
832                 for (zv = 0; zv < nVz; ++zv) {
833                   const PetscInt localVertex  = (zv*nVy +     0)*nVx + nVx-1 + nC; /* right bottom vertices */
834                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
835 
836                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
837                     localPoints[nL]        = localVertex;
838                     remotePoints[nL].rank  = neighbor;
839                     remotePoints[nL].index = remoteVertex;
840                     ++nL;
841                   }
842                 }
843               }
844             } else if (yp > 0) { /* top */
845               if (zp < 0) { /* back */
846                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
847                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
848 
849                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
850                   localPoints[nL]        = localVertex;
851                   remotePoints[nL].rank  = neighbor;
852                   remotePoints[nL].index = remoteVertex;
853                   ++nL;
854                 }
855               } else if (zp > 0) { /* front */
856                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
857                 const PetscInt remoteVertex = (      0*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
858 
859                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
860                   localPoints[nL]        = localVertex;
861                   remotePoints[nL].rank  = neighbor;
862                   remotePoints[nL].index = remoteVertex;
863                   ++nL;
864                 }
865               } else {
866                 for (zv = 0; zv < nVz; ++zv) {
867                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
868                   const PetscInt remoteVertex = (zv*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
869 
870                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
871                     localPoints[nL]        = localVertex;
872                     remotePoints[nL].rank  = neighbor;
873                     remotePoints[nL].index = remoteVertex;
874                     ++nL;
875                   }
876                 }
877               }
878             } else {
879               if (zp < 0) { /* back */
880                 for (yv = 0; yv < nVy; ++yv) {
881                   const PetscInt localVertex  = (      0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
882                   const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
883 
884                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
885                     localPoints[nL]        = localVertex;
886                     remotePoints[nL].rank  = neighbor;
887                     remotePoints[nL].index = remoteVertex;
888                     ++nL;
889                   }
890                 }
891               } else if (zp > 0) { /* front */
892                 for (yv = 0; yv < nVy; ++yv) {
893                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
894                   const PetscInt remoteVertex = (      0*nVy + yv)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */
895 
896                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
897                     localPoints[nL]        = localVertex;
898                     remotePoints[nL].rank  = neighbor;
899                     remotePoints[nL].index = remoteVertex;
900                     ++nL;
901                   }
902                 }
903               } else {
904                 for (zv = 0; zv < nVz; ++zv) {
905                   for (yv = 0; yv < nVy; ++yv) {
906                     const PetscInt localVertex  = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
907                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0     + nC; /* TODO: Correct this for neighbor sizes */
908 
909                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
910                       localPoints[nL]        = localVertex;
911                       remotePoints[nL].rank  = neighbor;
912                       remotePoints[nL].index = remoteVertex;
913                       ++nL;
914                     }
915                   }
916                 }
917 #if 0
918                 for (xf = 0; xf < nxF; ++xf) {
919                   /* THIS IS WRONG */
920                   const PetscInt localFace  = 0 + nC+nV; /* right faces */
921                   const PetscInt remoteFace = 0 + nC+nV;
922 
923                   if (!PetscBTLookupSet(isLeaf, localFace)) {
924                     localPoints[nL]        = localFace;
925                     remotePoints[nL].rank  = neighbor;
926                     remotePoints[nL].index = remoteFace;
927                     ++nL;
928                   }
929                 }
930 #endif
931               }
932             }
933           } else {
934             if (yp < 0) { /* bottom */
935               if (zp < 0) { /* back */
936                 for (xv = 0; xv < nVx; ++xv) {
937                   const PetscInt localVertex  = (      0*nVy +     0)*nVx + xv + nC; /* bottom back vertices */
938                   const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
939 
940                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
941                     localPoints[nL]        = localVertex;
942                     remotePoints[nL].rank  = neighbor;
943                     remotePoints[nL].index = remoteVertex;
944                     ++nL;
945                   }
946                 }
947               } else if (zp > 0) { /* front */
948                 for (xv = 0; xv < nVx; ++xv) {
949                   const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* bottom front vertices */
950                   const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
951 
952                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
953                     localPoints[nL]        = localVertex;
954                     remotePoints[nL].rank  = neighbor;
955                     remotePoints[nL].index = remoteVertex;
956                     ++nL;
957                   }
958                 }
959               } else {
960                 for (zv = 0; zv < nVz; ++zv) {
961                   for (xv = 0; xv < nVx; ++xv) {
962                     const PetscInt localVertex  = (zv*nVy +     0)*nVx + xv + nC; /* bottom vertices */
963                     const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
964 
965                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
966                       localPoints[nL]        = localVertex;
967                       remotePoints[nL].rank  = neighbor;
968                       remotePoints[nL].index = remoteVertex;
969                       ++nL;
970                     }
971                   }
972                 }
973 #if 0
974                 for (yf = 0; yf < nyF; ++yf) {
975                   /* THIS IS WRONG */
976                   const PetscInt localFace  = 0 + nC+nV; /* bottom faces */
977                   const PetscInt remoteFace = 0 + nC+nV;
978 
979                   if (!PetscBTLookupSet(isLeaf, localFace)) {
980                     localPoints[nL]        = localFace;
981                     remotePoints[nL].rank  = neighbor;
982                     remotePoints[nL].index = remoteFace;
983                     ++nL;
984                   }
985                 }
986 #endif
987               }
988             } else if (yp > 0) { /* top */
989               if (zp < 0) { /* back */
990                 for (xv = 0; xv < nVx; ++xv) {
991                   const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
992                   const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
993 
994                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
995                     localPoints[nL]        = localVertex;
996                     remotePoints[nL].rank  = neighbor;
997                     remotePoints[nL].index = remoteVertex;
998                     ++nL;
999                   }
1000                 }
1001               } else if (zp > 0) { /* front */
1002                 for (xv = 0; xv < nVx; ++xv) {
1003                   const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
1004                   const PetscInt remoteVertex = (      0*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1005 
1006                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
1007                     localPoints[nL]        = localVertex;
1008                     remotePoints[nL].rank  = neighbor;
1009                     remotePoints[nL].index = remoteVertex;
1010                     ++nL;
1011                   }
1012                 }
1013               } else {
1014                 for (zv = 0; zv < nVz; ++zv) {
1015                   for (xv = 0; xv < nVx; ++xv) {
1016                     const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
1017                     const PetscInt remoteVertex = (zv*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1018 
1019                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
1020                       localPoints[nL]        = localVertex;
1021                       remotePoints[nL].rank  = neighbor;
1022                       remotePoints[nL].index = remoteVertex;
1023                       ++nL;
1024                     }
1025                   }
1026                 }
1027 #if 0
1028                 for (yf = 0; yf < nyF; ++yf) {
1029                   /* THIS IS WRONG */
1030                   const PetscInt localFace  = 0 + nC+nV; /* top faces */
1031                   const PetscInt remoteFace = 0 + nC+nV;
1032 
1033                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1034                     localPoints[nL]        = localFace;
1035                     remotePoints[nL].rank  = neighbor;
1036                     remotePoints[nL].index = remoteFace;
1037                     ++nL;
1038                   }
1039                 }
1040 #endif
1041               }
1042             } else {
1043               if (zp < 0) { /* back */
1044                 for (yv = 0; yv < nVy; ++yv) {
1045                   for (xv = 0; xv < nVx; ++xv) {
1046                     const PetscInt localVertex  = (      0*nVy + yv)*nVx + xv + nC; /* back vertices */
1047                     const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1048 
1049                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
1050                       localPoints[nL]        = localVertex;
1051                       remotePoints[nL].rank  = neighbor;
1052                       remotePoints[nL].index = remoteVertex;
1053                       ++nL;
1054                     }
1055                   }
1056                 }
1057 #if 0
1058                 for (zf = 0; zf < nzF; ++zf) {
1059                   /* THIS IS WRONG */
1060                   const PetscInt localFace  = 0 + nC+nV; /* back faces */
1061                   const PetscInt remoteFace = 0 + nC+nV;
1062 
1063                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1064                     localPoints[nL]        = localFace;
1065                     remotePoints[nL].rank  = neighbor;
1066                     remotePoints[nL].index = remoteFace;
1067                     ++nL;
1068                   }
1069                 }
1070 #endif
1071               } else if (zp > 0) { /* front */
1072                 for (yv = 0; yv < nVy; ++yv) {
1073                   for (xv = 0; xv < nVx; ++xv) {
1074                     const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
1075                     const PetscInt remoteVertex = (      0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1076 
1077                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
1078                       localPoints[nL]        = localVertex;
1079                       remotePoints[nL].rank  = neighbor;
1080                       remotePoints[nL].index = remoteVertex;
1081                       ++nL;
1082                     }
1083                   }
1084                 }
1085 #if 0
1086                 for (zf = 0; zf < nzF; ++zf) {
1087                   /* THIS IS WRONG */
1088                   const PetscInt localFace  = 0 + nC+nV; /* front faces */
1089                   const PetscInt remoteFace = 0 + nC+nV;
1090 
1091                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1092                     localPoints[nL]        = localFace;
1093                     remotePoints[nL].rank  = neighbor;
1094                     remotePoints[nL].index = remoteFace;
1095                     ++nL;
1096                   }
1097                 }
1098 #endif
1099               } else {
1100                 /* Nothing is shared from the interior */
1101               }
1102             }
1103           }
1104         }
1105       }
1106     }
1107   }
1108   ierr = PetscBTDestroy(&isLeaf);CHKERRQ(ierr);
1109   /* Remove duplication in leaf determination */
1110   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);
1111   ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);CHKERRQ(ierr);
1112   ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr);
1113   ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr);
1114   ierr = PetscSFDestroy(&sf);CHKERRQ(ierr);
1115   ierr = DMSetDefaultSection(dm, section);CHKERRQ(ierr);
1116 #undef __FUNCT__
1117 #define __FUNCT__ "DMDASetVertexCoordinates"
1118 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
1119 {
1120   DM_DA         *da = (DM_DA *) dm->data;
1121   Vec            coordinates;
1122   PetscSection   section;
1123   PetscScalar   *coords;
1124   PetscReal      h[3];
1125   PetscInt       dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
1126   PetscErrorCode ierr;
1127 
1128   PetscFunctionBegin;
1129   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1130   ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
1131   h[0] = (xu - xl)/M;
1132   h[1] = (yu - yl)/N;
1133   h[2] = (zu - zl)/P;
1134   ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
1135   ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr);
1136   ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);CHKERRQ(ierr);
1137   ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr);
1138   ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr);
1139   ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr);
1140   for (v = vStart; v < vEnd; ++v) {
1141     ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr);
1142   }
1143   ierr = PetscSectionSetUp(section);CHKERRQ(ierr);
1144   ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
1145   ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr);
1146   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
1147   for (k = 0; k < nVz; ++k) {
1148     PetscInt ind[3] = {0, 0, k + da->zs}, d, off;
1149 
1150     for (j = 0; j < nVy; ++j) {
1151       ind[1] = j + da->ys;
1152       for (i = 0; i < nVx; ++i) {
1153         const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;
1154 
1155         ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr);
1156         ind[0] = i + da->xs;
1157         for (d = 0; d < dim; ++d) {
1158           coords[off+d] = h[d]*ind[d];
1159         }
1160       }
1161     }
1162   }
1163   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
1164   ierr = DMSetCoordinateSection(dm, section);CHKERRQ(ierr);
1165   ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
1166   ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);
1167   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
1168   PetscFunctionReturn(0);
1169 }
1170   PetscFunctionReturn(0);
1171 }
1172 
1173 /* ------------------------------------------------------------------- */
1174 
1175 #undef __FUNCT__
1176 #define __FUNCT__ "DMDAGetArray"
1177 /*@C
1178      DMDAGetArray - Gets a work array for a DMDA
1179 
1180     Input Parameter:
1181 +    da - information about my local patch
1182 -    ghosted - do you want arrays for the ghosted or nonghosted patch
1183 
1184     Output Parameters:
1185 .    vptr - array data structured
1186 
1187     Note:  The vector values are NOT initialized and may have garbage in them, so you may need
1188            to zero them.
1189 
1190   Level: advanced
1191 
1192 .seealso: DMDARestoreArray()
1193 
1194 @*/
1195 PetscErrorCode  DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
1196 {
1197   PetscErrorCode ierr;
1198   PetscInt       j,i,xs,ys,xm,ym,zs,zm;
1199   char           *iarray_start;
1200   void           **iptr = (void**)vptr;
1201   DM_DA          *dd    = (DM_DA*)da->data;
1202 
1203   PetscFunctionBegin;
1204   PetscValidHeaderSpecific(da,DM_CLASSID,1);
1205   if (ghosted) {
1206     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1207       if (dd->arrayghostedin[i]) {
1208         *iptr                 = dd->arrayghostedin[i];
1209         iarray_start          = (char*)dd->startghostedin[i];
1210         dd->arrayghostedin[i] = NULL;
1211         dd->startghostedin[i] = NULL;
1212 
1213         goto done;
1214       }
1215     }
1216     xs = dd->Xs;
1217     ys = dd->Ys;
1218     zs = dd->Zs;
1219     xm = dd->Xe-dd->Xs;
1220     ym = dd->Ye-dd->Ys;
1221     zm = dd->Ze-dd->Zs;
1222   } else {
1223     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1224       if (dd->arrayin[i]) {
1225         *iptr          = dd->arrayin[i];
1226         iarray_start   = (char*)dd->startin[i];
1227         dd->arrayin[i] = NULL;
1228         dd->startin[i] = NULL;
1229 
1230         goto done;
1231       }
1232     }
1233     xs = dd->xs;
1234     ys = dd->ys;
1235     zs = dd->zs;
1236     xm = dd->xe-dd->xs;
1237     ym = dd->ye-dd->ys;
1238     zm = dd->ze-dd->zs;
1239   }
1240 
1241   switch (dd->dim) {
1242   case 1: {
1243     void *ptr;
1244 
1245     ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
1246 
1247     ptr   = (void*)(iarray_start - xs*sizeof(PetscScalar));
1248     *iptr = (void*)ptr;
1249     break;
1250   }
1251   case 2: {
1252     void **ptr;
1253 
1254     ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
1255 
1256     ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
1257     for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
1258     *iptr = (void*)ptr;
1259     break;
1260   }
1261   case 3: {
1262     void ***ptr,**bptr;
1263 
1264     ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr);
1265 
1266     ptr  = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1267     bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1268     for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1269     for (i=zs; i<zs+zm; i++) {
1270       for (j=ys; j<ys+ym; j++) {
1271         ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1272       }
1273     }
1274     *iptr = (void*)ptr;
1275     break;
1276   }
1277   default:
1278     SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
1279   }
1280 
1281 done:
1282   /* add arrays to the checked out list */
1283   if (ghosted) {
1284     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1285       if (!dd->arrayghostedout[i]) {
1286         dd->arrayghostedout[i] = *iptr;
1287         dd->startghostedout[i] = iarray_start;
1288         break;
1289       }
1290     }
1291   } else {
1292     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1293       if (!dd->arrayout[i]) {
1294         dd->arrayout[i] = *iptr;
1295         dd->startout[i] = iarray_start;
1296         break;
1297       }
1298     }
1299   }
1300   PetscFunctionReturn(0);
1301 }
1302 
1303 #undef __FUNCT__
1304 #define __FUNCT__ "DMDARestoreArray"
1305 /*@C
1306      DMDARestoreArray - Restores an array of derivative types for a DMDA
1307 
1308     Input Parameter:
1309 +    da - information about my local patch
1310 .    ghosted - do you want arrays for the ghosted or nonghosted patch
1311 -    vptr - array data structured to be passed to ad_FormFunctionLocal()
1312 
1313      Level: advanced
1314 
1315 .seealso: DMDAGetArray()
1316 
1317 @*/
1318 PetscErrorCode  DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1319 {
1320   PetscInt i;
1321   void     **iptr = (void**)vptr,*iarray_start = 0;
1322   DM_DA    *dd    = (DM_DA*)da->data;
1323 
1324   PetscFunctionBegin;
1325   PetscValidHeaderSpecific(da,DM_CLASSID,1);
1326   if (ghosted) {
1327     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1328       if (dd->arrayghostedout[i] == *iptr) {
1329         iarray_start           = dd->startghostedout[i];
1330         dd->arrayghostedout[i] = NULL;
1331         dd->startghostedout[i] = NULL;
1332         break;
1333       }
1334     }
1335     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1336       if (!dd->arrayghostedin[i]) {
1337         dd->arrayghostedin[i] = *iptr;
1338         dd->startghostedin[i] = iarray_start;
1339         break;
1340       }
1341     }
1342   } else {
1343     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1344       if (dd->arrayout[i] == *iptr) {
1345         iarray_start    = dd->startout[i];
1346         dd->arrayout[i] = NULL;
1347         dd->startout[i] = NULL;
1348         break;
1349       }
1350     }
1351     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1352       if (!dd->arrayin[i]) {
1353         dd->arrayin[i] = *iptr;
1354         dd->startin[i] = iarray_start;
1355         break;
1356       }
1357     }
1358   }
1359   PetscFunctionReturn(0);
1360 }
1361 
1362