Actual source code: dalocal.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  2: /*
  3:   Code for manipulating distributed regular arrays in parallel.
  4: */

  6:  #include <petsc/private/dmdaimpl.h>
  7:  #include <petscbt.h>
  8:  #include <petscsf.h>
  9:  #include <petscds.h>
 10:  #include <petscfe.h>

 12: /*
 13:    This allows the DMDA vectors to properly tell MATLAB their dimensions
 14: */
 15: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 16: #include <engine.h>   /* MATLAB include file */
 17: #include <mex.h>      /* MATLAB include file */
 18: static PetscErrorCode  VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
 19: {
 21:   PetscInt       n,m;
 22:   Vec            vec = (Vec)obj;
 23:   PetscScalar    *array;
 24:   mxArray        *mat;
 25:   DM             da;

 28:   VecGetDM(vec, &da);
 29:   if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA");
 30:   DMDAGetGhostCorners(da,0,0,0,&m,&n,0);

 32:   VecGetArray(vec,&array);
 33: #if !defined(PETSC_USE_COMPLEX)
 34:   mat = mxCreateDoubleMatrix(m,n,mxREAL);
 35: #else
 36:   mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
 37: #endif
 38:   PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));
 39:   PetscObjectName(obj);
 40:   engPutVariable((Engine*)mengine,obj->name,mat);

 42:   VecRestoreArray(vec,&array);
 43:   return(0);
 44: }
 45: #endif


 48: PetscErrorCode  DMCreateLocalVector_DA(DM da,Vec *g)
 49: {
 51:   DM_DA          *dd = (DM_DA*)da->data;

 56:   if (da->defaultSection) {
 57:     DMCreateLocalVector_Section_Private(da,g);
 58:   } else {
 59:     VecCreate(PETSC_COMM_SELF,g);
 60:     VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);
 61:     VecSetBlockSize(*g,dd->w);
 62:     VecSetType(*g,da->vectype);
 63:     VecSetDM(*g, da);
 64: #if defined(PETSC_HAVE_MATLAB_ENGINE)
 65:     if (dd->w == 1  && da->dim == 2) {
 66:       PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);
 67:     }
 68: #endif
 69:   }
 70:   return(0);
 71: }

 73: /*@
 74:   DMDAGetNumCells - Get the number of cells in the local piece of the DMDA. This includes ghost cells.

 76:   Input Parameter:
 77: . dm - The DM object

 79:   Output Parameters:
 80: + numCellsX - The number of local cells in the x-direction
 81: . numCellsY - The number of local cells in the y-direction
 82: . numCellsZ - The number of local cells in the z-direction
 83: - numCells - The number of local cells

 85:   Level: developer

 87: .seealso: DMDAGetCellPoint()
 88: @*/
 89: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
 90: {
 91:   DM_DA         *da  = (DM_DA*) dm->data;
 92:   const PetscInt dim = dm->dim;
 93:   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
 94:   const PetscInt nC  = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);

 98:   if (numCellsX) {
100:     *numCellsX = mx;
101:   }
102:   if (numCellsY) {
104:     *numCellsY = my;
105:   }
106:   if (numCellsZ) {
108:     *numCellsZ = mz;
109:   }
110:   if (numCells) {
112:     *numCells = nC;
113:   }
114:   return(0);
115: }

117: /*@
118:   DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the DMDA

120:   Input Parameters:
121: + dm - The DM object
122: - i,j,k - The global indices for the cell

124:   Output Parameters:
125: . point - The local DM point

127:   Level: developer

129: .seealso: DMDAGetNumCells()
130: @*/
131: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
132: {
133:   const PetscInt dim = dm->dim;
134:   DMDALocalInfo  info;

140:   DMDAGetLocalInfo(dm, &info);
141:   if (dim > 0) {if ((i < info.gxs) || (i >= info.gxs+info.gxm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "X index %d not in [%d, %d)", i, info.gxs, info.gxs+info.gxm);}
142:   if (dim > 1) {if ((j < info.gys) || (j >= info.gys+info.gym)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %d not in [%d, %d)", j, info.gys, info.gys+info.gym);}
143:   if (dim > 2) {if ((k < info.gzs) || (k >= info.gzs+info.gzm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %d not in [%d, %d)", k, info.gzs, info.gzs+info.gzm);}
144:   *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0);
145:   return(0);
146: }

148: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
149: {
150:   DM_DA          *da = (DM_DA*) dm->data;
151:   const PetscInt dim = dm->dim;
152:   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
153:   const PetscInt nVx = mx+1;
154:   const PetscInt nVy = dim > 1 ? (my+1) : 1;
155:   const PetscInt nVz = dim > 2 ? (mz+1) : 1;
156:   const PetscInt nV  = nVx*nVy*nVz;

159:   if (numVerticesX) {
161:     *numVerticesX = nVx;
162:   }
163:   if (numVerticesY) {
165:     *numVerticesY = nVy;
166:   }
167:   if (numVerticesZ) {
169:     *numVerticesZ = nVz;
170:   }
171:   if (numVertices) {
173:     *numVertices = nV;
174:   }
175:   return(0);
176: }

178: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
179: {
180:   DM_DA          *da = (DM_DA*) dm->data;
181:   const PetscInt dim = dm->dim;
182:   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
183:   const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
184:   const PetscInt nXF = (mx+1)*nxF;
185:   const PetscInt nyF = mx*(dim > 2 ? mz : 1);
186:   const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0;
187:   const PetscInt nzF = mx*(dim > 1 ? my : 0);
188:   const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0;

191:   if (numXFacesX) {
193:     *numXFacesX = nxF;
194:   }
195:   if (numXFaces) {
197:     *numXFaces = nXF;
198:   }
199:   if (numYFacesY) {
201:     *numYFacesY = nyF;
202:   }
203:   if (numYFaces) {
205:     *numYFaces = nYF;
206:   }
207:   if (numZFacesZ) {
209:     *numZFacesZ = nzF;
210:   }
211:   if (numZFaces) {
213:     *numZFaces = nZF;
214:   }
215:   return(0);
216: }

218: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
219: {
220:   const PetscInt dim = dm->dim;
221:   PetscInt       nC, nV, nXF, nYF, nZF;

227:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
228:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
229:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
230:   if (height == 0) {
231:     /* Cells */
232:     if (pStart) *pStart = 0;
233:     if (pEnd)   *pEnd   = nC;
234:   } else if (height == 1) {
235:     /* Faces */
236:     if (pStart) *pStart = nC+nV;
237:     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
238:   } else if (height == dim) {
239:     /* Vertices */
240:     if (pStart) *pStart = nC;
241:     if (pEnd)   *pEnd   = nC+nV;
242:   } else if (height < 0) {
243:     /* All points */
244:     if (pStart) *pStart = 0;
245:     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
246:   } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
247:   return(0);
248: }

250: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
251: {
252:   const PetscInt dim = dm->dim;
253:   PetscInt       nC, nV, nXF, nYF, nZF;

259:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
260:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
261:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
262:   if (depth == dim) {
263:     /* Cells */
264:     if (pStart) *pStart = 0;
265:     if (pEnd)   *pEnd   = nC;
266:   } else if (depth == dim-1) {
267:     /* Faces */
268:     if (pStart) *pStart = nC+nV;
269:     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
270:   } else if (depth == 0) {
271:     /* Vertices */
272:     if (pStart) *pStart = nC;
273:     if (pEnd)   *pEnd   = nC+nV;
274:   } else if (depth < 0) {
275:     /* All points */
276:     if (pStart) *pStart = 0;
277:     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
278:   } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
279:   return(0);
280: }

282: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
283: {
284:   const PetscInt dim = dm->dim;
285:   PetscInt       nC, nV, nXF, nYF, nZF;

289:   *coneSize = 0;
290:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
291:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
292:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
293:   switch (dim) {
294:   case 2:
295:     if (p >= 0) {
296:       if (p < nC) {
297:         *coneSize = 4;
298:       } else if (p < nC+nV) {
299:         *coneSize = 0;
300:       } else if (p < nC+nV+nXF+nYF+nZF) {
301:         *coneSize = 2;
302:       } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
303:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
304:     break;
305:   case 3:
306:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
307:     break;
308:   }
309:   return(0);
310: }

312: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
313: {
314:   const PetscInt dim = dm->dim;
315:   PetscInt       nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;

319:   if (!cone) {DMGetWorkArray(dm, 6, MPIU_INT, cone);}
320:   DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);
321:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
322:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
323:   switch (dim) {
324:   case 2:
325:     if (p >= 0) {
326:       if (p < nC) {
327:         const PetscInt cy = p / nCx;
328:         const PetscInt cx = p % nCx;

330:         (*cone)[0] = cy*nxF + cx + nC+nV;
331:         (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
332:         (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
333:         (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
334:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
335:       } else if (p < nC+nV) {
336:       } else if (p < nC+nV+nXF) {
337:         const PetscInt fy = (p - nC+nV) / nxF;
338:         const PetscInt fx = (p - nC+nV) % nxF;

340:         (*cone)[0] = fy*nVx + fx + nC;
341:         (*cone)[1] = fy*nVx + fx + 1 + nC;
342:       } else if (p < nC+nV+nXF+nYF) {
343:         const PetscInt fx = (p - nC+nV+nXF) / nyF;
344:         const PetscInt fy = (p - nC+nV+nXF) % nyF;

346:         (*cone)[0] = fy*nVx + fx + nC;
347:         (*cone)[1] = fy*nVx + fx + nVx + nC;
348:       } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
349:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
350:     break;
351:   case 3:
352:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
353:     break;
354:   }
355:   return(0);
356: }

358: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
359: {

363:   DMGetWorkArray(dm, 6, MPIU_INT, cone);
364:   return(0);
365: }

367: /*@C
368:   DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with
369:   different numbers of dofs on vertices, cells, and faces in each direction.

371:   Input Parameters:
372: + dm- The DMDA
373: . numFields - The number of fields
374: . numComp - The number of components in each field
375: . numDof - The number of dofs per dimension for each field
376: . numFaceDof - The number of dofs per face for each field and direction, or NULL

378:   Level: developer

380:   Note:
381:   The default DMDA numbering is as follows:

383:     - Cells:    [0,             nC)
384:     - Vertices: [nC,            nC+nV)
385:     - X-Faces:  [nC+nV,         nC+nV+nXF)         normal is +- x-dir
386:     - Y-Faces:  [nC+nV+nXF,     nC+nV+nXF+nYF)     normal is +- y-dir
387:     - Z-Faces:  [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir

389:   We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment.
390: @*/
391: PetscErrorCode DMDACreateSection(DM dm, const PetscInt numComp[], const PetscInt numDof[], const PetscInt numFaceDof[], PetscSection *s)
392: {
393:   PetscSection      section;
394:   const PetscInt    dim = dm->dim;
395:   PetscInt          numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0};
396:   PetscBT           isLeaf;
397:   PetscSF           sf;
398:   PetscMPIInt       rank;
399:   const PetscMPIInt *neighbors;
400:   PetscInt          *localPoints;
401:   PetscSFNode       *remotePoints;
402:   PetscInt          nleaves = 0,  nleavesCheck = 0, nL = 0;
403:   PetscInt          nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF;
404:   PetscInt          pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd;
405:   PetscInt          f, v, c, xf, yf, zf, xn, yn, zn;
406:   PetscErrorCode    ierr;

411:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
412:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
413:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
414:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
415:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
416:   DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);
417:   DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);
418:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
419:   xfStart = vEnd;  xfEnd = xfStart+nXF;
420:   yfStart = xfEnd; yfEnd = yfStart+nYF;
421:   zfStart = yfEnd; zfEnd = zfStart+nZF;
422:   if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd);
423:   /* Create local section */
424:   DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);
425:   for (f = 0; f < numFields; ++f) {
426:     numVertexTotDof  += numDof[f*(dim+1)+0];
427:     numCellTotDof    += numDof[f*(dim+1)+dim];
428:     numFaceTotDof[0] += dim > 0 ? (numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]) : 0;
429:     numFaceTotDof[1] += dim > 1 ? (numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]) : 0;
430:     numFaceTotDof[2] += dim > 2 ? (numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]) : 0;
431:   }
432:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
433:   if (numFields > 0) {
434:     PetscSectionSetNumFields(section, numFields);
435:     for (f = 0; f < numFields; ++f) {
436:       const char *name;

438:       DMDAGetFieldName(dm, f, &name);
439:       PetscSectionSetFieldName(section, f, name ? name : "Field");
440:       if (numComp) {
441:         PetscSectionSetFieldComponents(section, f, numComp[f]);
442:       }
443:     }
444:   }
445:   PetscSectionSetChart(section, pStart, pEnd);
446:   for (v = vStart; v < vEnd; ++v) {
447:     for (f = 0; f < numFields; ++f) {
448:       PetscSectionSetFieldDof(section, v, f, numDof[f*(dim+1)+0]);
449:     }
450:     PetscSectionSetDof(section, v, numVertexTotDof);
451:   }
452:   for (xf = xfStart; xf < xfEnd; ++xf) {
453:     for (f = 0; f < numFields; ++f) {
454:       PetscSectionSetFieldDof(section, xf, f, numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]);
455:     }
456:     PetscSectionSetDof(section, xf, numFaceTotDof[0]);
457:   }
458:   for (yf = yfStart; yf < yfEnd; ++yf) {
459:     for (f = 0; f < numFields; ++f) {
460:       PetscSectionSetFieldDof(section, yf, f, numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]);
461:     }
462:     PetscSectionSetDof(section, yf, numFaceTotDof[1]);
463:   }
464:   for (zf = zfStart; zf < zfEnd; ++zf) {
465:     for (f = 0; f < numFields; ++f) {
466:       PetscSectionSetFieldDof(section, zf, f, numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]);
467:     }
468:     PetscSectionSetDof(section, zf, numFaceTotDof[2]);
469:   }
470:   for (c = cStart; c < cEnd; ++c) {
471:     for (f = 0; f < numFields; ++f) {
472:       PetscSectionSetFieldDof(section, c, f, numDof[f*(dim+1)+dim]);
473:     }
474:     PetscSectionSetDof(section, c, numCellTotDof);
475:   }
476:   PetscSectionSetUp(section);
477:   /* Create mesh point SF */
478:   PetscBTCreate(pEnd-pStart, &isLeaf);
479:   DMDAGetNeighbors(dm, &neighbors);
480:   for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
481:     for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
482:       for (xn = 0; xn < 3; ++xn) {
483:         const PetscInt xp       = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
484:         const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
485:         PetscInt       xv, yv, zv;

487:         if (neighbor >= 0 && neighbor < rank) {
488:           if (xp < 0) { /* left */
489:             if (yp < 0) { /* bottom */
490:               if (zp < 0) { /* back */
491:                 const PetscInt localVertex  = (      0*nVy +     0)*nVx +     0 + nC; /* left bottom back vertex */
492:                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
493:               } else if (zp > 0) { /* front */
494:                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* left bottom front vertex */
495:                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
496:               } else {
497:                 for (zv = 0; zv < nVz; ++zv) {
498:                   const PetscInt localVertex  = (zv*nVy +     0)*nVx +     0 + nC; /* left bottom vertices */
499:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
500:                 }
501:               }
502:             } else if (yp > 0) { /* top */
503:               if (zp < 0) { /* back */
504:                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx +     0 + nC; /* left top back vertex */
505:                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
506:               } else if (zp > 0) { /* front */
507:                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* left top front vertex */
508:                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
509:               } else {
510:                 for (zv = 0; zv < nVz; ++zv) {
511:                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx +     0 + nC; /* left top vertices */
512:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
513:                 }
514:               }
515:             } else {
516:               if (zp < 0) { /* back */
517:                 for (yv = 0; yv < nVy; ++yv) {
518:                   const PetscInt localVertex  = (      0*nVy + yv)*nVx +     0 + nC; /* left back vertices */
519:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
520:                 }
521:               } else if (zp > 0) { /* front */
522:                 for (yv = 0; yv < nVy; ++yv) {
523:                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx +     0 + nC; /* left front vertices */
524:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
525:                 }
526:               } else {
527:                 for (zv = 0; zv < nVz; ++zv) {
528:                   for (yv = 0; yv < nVy; ++yv) {
529:                     const PetscInt localVertex  = (zv*nVy + yv)*nVx +     0 + nC; /* left vertices */
530:                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
531:                   }
532:                 }
533: #if 0
534:                 for (xf = 0; xf < nxF; ++xf) {
535:                   /* THIS IS WRONG */
536:                   const PetscInt localFace  = 0 + nC+nV; /* left faces */
537:                   if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
538:                 }
539: #endif
540:               }
541:             }
542:           } else if (xp > 0) { /* right */
543:             if (yp < 0) { /* bottom */
544:               if (zp < 0) { /* back */
545:                 const PetscInt localVertex  = (      0*nVy +     0)*nVx + nVx-1 + nC; /* right bottom back vertex */
546:                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
547:               } else if (zp > 0) { /* front */
548:                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* right bottom front vertex */
549:                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
550:               } else {
551:                 for (zv = 0; zv < nVz; ++zv) {
552:                   const PetscInt localVertex  = (zv*nVy +     0)*nVx + nVx-1 + nC; /* right bottom vertices */
553:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
554:                 }
555:               }
556:             } else if (yp > 0) { /* top */
557:               if (zp < 0) { /* back */
558:                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
559:                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
560:               } else if (zp > 0) { /* front */
561:                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
562:                 if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
563:               } else {
564:                 for (zv = 0; zv < nVz; ++zv) {
565:                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
566:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
567:                 }
568:               }
569:             } else {
570:               if (zp < 0) { /* back */
571:                 for (yv = 0; yv < nVy; ++yv) {
572:                   const PetscInt localVertex  = (      0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
573:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
574:                 }
575:               } else if (zp > 0) { /* front */
576:                 for (yv = 0; yv < nVy; ++yv) {
577:                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
578:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
579:                 }
580:               } else {
581:                 for (zv = 0; zv < nVz; ++zv) {
582:                   for (yv = 0; yv < nVy; ++yv) {
583:                     const PetscInt localVertex  = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
584:                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
585:                   }
586:                 }
587: #if 0
588:                 for (xf = 0; xf < nxF; ++xf) {
589:                   /* THIS IS WRONG */
590:                   const PetscInt localFace  = 0 + nC+nV; /* right faces */
591:                   if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
592:                 }
593: #endif
594:               }
595:             }
596:           } else {
597:             if (yp < 0) { /* bottom */
598:               if (zp < 0) { /* back */
599:                 for (xv = 0; xv < nVx; ++xv) {
600:                   const PetscInt localVertex  = (      0*nVy +     0)*nVx + xv + nC; /* bottom back vertices */
601:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
602:                 }
603:               } else if (zp > 0) { /* front */
604:                 for (xv = 0; xv < nVx; ++xv) {
605:                   const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* bottom front vertices */
606:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
607:                 }
608:               } else {
609:                 for (zv = 0; zv < nVz; ++zv) {
610:                   for (xv = 0; xv < nVx; ++xv) {
611:                     const PetscInt localVertex  = (zv*nVy +     0)*nVx + xv + nC; /* bottom vertices */
612:                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
613:                   }
614:                 }
615: #if 0
616:                 for (yf = 0; yf < nyF; ++yf) {
617:                   /* THIS IS WRONG */
618:                   const PetscInt localFace  = 0 + nC+nV; /* bottom faces */
619:                   if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves;
620:                 }
621: #endif
622:               }
623:             } else if (yp > 0) { /* top */
624:               if (zp < 0) { /* back */
625:                 for (xv = 0; xv < nVx; ++xv) {
626:                   const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
627:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
628:                 }
629:               } else if (zp > 0) { /* front */
630:                 for (xv = 0; xv < nVx; ++xv) {
631:                   const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
632:                   if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
633:                 }
634:               } else {
635:                 for (zv = 0; zv < nVz; ++zv) {
636:                   for (xv = 0; xv < nVx; ++xv) {
637:                     const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
638:                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
639:                   }
640:                 }
641: #if 0
642:                 for (yf = 0; yf < nyF; ++yf) {
643:                   /* THIS IS WRONG */
644:                   const PetscInt localFace  = 0 + nC+nV; /* top faces */
645:                   if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves;
646:                 }
647: #endif
648:               }
649:             } else {
650:               if (zp < 0) { /* back */
651:                 for (yv = 0; yv < nVy; ++yv) {
652:                   for (xv = 0; xv < nVx; ++xv) {
653:                     const PetscInt localVertex  = (      0*nVy + yv)*nVx + xv + nC; /* back vertices */
654:                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
655:                   }
656:                 }
657: #if 0
658:                 for (zf = 0; zf < nzF; ++zf) {
659:                   /* THIS IS WRONG */
660:                   const PetscInt localFace  = 0 + nC+nV; /* back faces */
661:                   if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
662:                 }
663: #endif
664:               } else if (zp > 0) { /* front */
665:                 for (yv = 0; yv < nVy; ++yv) {
666:                   for (xv = 0; xv < nVx; ++xv) {
667:                     const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
668:                     if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
669:                   }
670:                 }
671: #if 0
672:                 for (zf = 0; zf < nzF; ++zf) {
673:                   /* THIS IS WRONG */
674:                   const PetscInt localFace  = 0 + nC+nV; /* front faces */
675:                   if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
676:                 }
677: #endif
678:               } else {
679:                 /* Nothing is shared from the interior */
680:               }
681:             }
682:           }
683:         }
684:       }
685:     }
686:   }
687:   PetscBTMemzero(pEnd-pStart, isLeaf);
688:   PetscMalloc1(nleaves,&localPoints);
689:   PetscMalloc1(nleaves,&remotePoints);
690:   for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
691:     for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
692:       for (xn = 0; xn < 3; ++xn) {
693:         const PetscInt xp       = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
694:         const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
695:         PetscInt       xv, yv, zv;

697:         if (neighbor >= 0 && neighbor < rank) {
698:           if (xp < 0) { /* left */
699:             if (yp < 0) { /* bottom */
700:               if (zp < 0) { /* back */
701:                 const PetscInt localVertex  = (      0*nVy +     0)*nVx +     0 + nC; /* left bottom back vertex */
702:                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

704:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
705:                   localPoints[nL]        = localVertex;
706:                   remotePoints[nL].rank  = neighbor;
707:                   remotePoints[nL].index = remoteVertex;
708:                   ++nL;
709:                 }
710:               } else if (zp > 0) { /* front */
711:                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* left bottom front vertex */
712:                 const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

714:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
715:                   localPoints[nL]        = localVertex;
716:                   remotePoints[nL].rank  = neighbor;
717:                   remotePoints[nL].index = remoteVertex;
718:                   ++nL;
719:                 }
720:               } else {
721:                 for (zv = 0; zv < nVz; ++zv) {
722:                   const PetscInt localVertex  = (zv*nVy +     0)*nVx +     0 + nC; /* left bottom vertices */
723:                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

725:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
726:                     localPoints[nL]        = localVertex;
727:                     remotePoints[nL].rank  = neighbor;
728:                     remotePoints[nL].index = remoteVertex;
729:                     ++nL;
730:                   }
731:                 }
732:               }
733:             } else if (yp > 0) { /* top */
734:               if (zp < 0) { /* back */
735:                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx +     0 + nC; /* left top back vertex */
736:                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

738:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
739:                   localPoints[nL]        = localVertex;
740:                   remotePoints[nL].rank  = neighbor;
741:                   remotePoints[nL].index = remoteVertex;
742:                   ++nL;
743:                 }
744:               } else if (zp > 0) { /* front */
745:                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* left top front vertex */
746:                 const PetscInt remoteVertex = (      0*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

748:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
749:                   localPoints[nL]        = localVertex;
750:                   remotePoints[nL].rank  = neighbor;
751:                   remotePoints[nL].index = remoteVertex;
752:                   ++nL;
753:                 }
754:               } else {
755:                 for (zv = 0; zv < nVz; ++zv) {
756:                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx +     0 + nC; /* left top vertices */
757:                   const PetscInt remoteVertex = (zv*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

759:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
760:                     localPoints[nL]        = localVertex;
761:                     remotePoints[nL].rank  = neighbor;
762:                     remotePoints[nL].index = remoteVertex;
763:                     ++nL;
764:                   }
765:                 }
766:               }
767:             } else {
768:               if (zp < 0) { /* back */
769:                 for (yv = 0; yv < nVy; ++yv) {
770:                   const PetscInt localVertex  = (      0*nVy + yv)*nVx +     0 + nC; /* left back vertices */
771:                   const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

773:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
774:                     localPoints[nL]        = localVertex;
775:                     remotePoints[nL].rank  = neighbor;
776:                     remotePoints[nL].index = remoteVertex;
777:                     ++nL;
778:                   }
779:                 }
780:               } else if (zp > 0) { /* front */
781:                 for (yv = 0; yv < nVy; ++yv) {
782:                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx +     0 + nC; /* left front vertices */
783:                   const PetscInt remoteVertex = (      0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

785:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
786:                     localPoints[nL]        = localVertex;
787:                     remotePoints[nL].rank  = neighbor;
788:                     remotePoints[nL].index = remoteVertex;
789:                     ++nL;
790:                   }
791:                 }
792:               } else {
793:                 for (zv = 0; zv < nVz; ++zv) {
794:                   for (yv = 0; yv < nVy; ++yv) {
795:                     const PetscInt localVertex  = (zv*nVy + yv)*nVx +     0 + nC; /* left vertices */
796:                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

798:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
799:                       localPoints[nL]        = localVertex;
800:                       remotePoints[nL].rank  = neighbor;
801:                       remotePoints[nL].index = remoteVertex;
802:                       ++nL;
803:                     }
804:                   }
805:                 }
806: #if 0
807:                 for (xf = 0; xf < nxF; ++xf) {
808:                   /* THIS IS WRONG */
809:                   const PetscInt localFace  = 0 + nC+nV; /* left faces */
810:                   const PetscInt remoteFace = 0 + nC+nV;

812:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
813:                     localPoints[nL]        = localFace;
814:                     remotePoints[nL].rank  = neighbor;
815:                     remotePoints[nL].index = remoteFace;
816:                   }
817:                 }
818: #endif
819:               }
820:             }
821:           } else if (xp > 0) { /* right */
822:             if (yp < 0) { /* bottom */
823:               if (zp < 0) { /* back */
824:                 const PetscInt localVertex  = (      0*nVy +     0)*nVx + nVx-1 + nC; /* right bottom back vertex */
825:                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

827:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
828:                   localPoints[nL]        = localVertex;
829:                   remotePoints[nL].rank  = neighbor;
830:                   remotePoints[nL].index = remoteVertex;
831:                   ++nL;
832:                 }
833:               } else if (zp > 0) { /* front */
834:                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* right bottom front vertex */
835:                 const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

837:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
838:                   localPoints[nL]        = localVertex;
839:                   remotePoints[nL].rank  = neighbor;
840:                   remotePoints[nL].index = remoteVertex;
841:                   ++nL;
842:                 }
843:               } else {
844:                 nleavesCheck += nVz;
845:                 for (zv = 0; zv < nVz; ++zv) {
846:                   const PetscInt localVertex  = (zv*nVy +     0)*nVx + nVx-1 + nC; /* right bottom vertices */
847:                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

849:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
850:                     localPoints[nL]        = localVertex;
851:                     remotePoints[nL].rank  = neighbor;
852:                     remotePoints[nL].index = remoteVertex;
853:                     ++nL;
854:                   }
855:                 }
856:               }
857:             } else if (yp > 0) { /* top */
858:               if (zp < 0) { /* back */
859:                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
860:                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

862:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
863:                   localPoints[nL]        = localVertex;
864:                   remotePoints[nL].rank  = neighbor;
865:                   remotePoints[nL].index = remoteVertex;
866:                   ++nL;
867:                 }
868:               } else if (zp > 0) { /* front */
869:                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
870:                 const PetscInt remoteVertex = (      0*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

872:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
873:                   localPoints[nL]        = localVertex;
874:                   remotePoints[nL].rank  = neighbor;
875:                   remotePoints[nL].index = remoteVertex;
876:                   ++nL;
877:                 }
878:               } else {
879:                 for (zv = 0; zv < nVz; ++zv) {
880:                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
881:                   const PetscInt remoteVertex = (zv*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

883:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
884:                     localPoints[nL]        = localVertex;
885:                     remotePoints[nL].rank  = neighbor;
886:                     remotePoints[nL].index = remoteVertex;
887:                     ++nL;
888:                   }
889:                 }
890:               }
891:             } else {
892:               if (zp < 0) { /* back */
893:                 for (yv = 0; yv < nVy; ++yv) {
894:                   const PetscInt localVertex  = (      0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
895:                   const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

897:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
898:                     localPoints[nL]        = localVertex;
899:                     remotePoints[nL].rank  = neighbor;
900:                     remotePoints[nL].index = remoteVertex;
901:                     ++nL;
902:                   }
903:                 }
904:               } else if (zp > 0) { /* front */
905:                 for (yv = 0; yv < nVy; ++yv) {
906:                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
907:                   const PetscInt remoteVertex = (      0*nVy + yv)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

909:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
910:                     localPoints[nL]        = localVertex;
911:                     remotePoints[nL].rank  = neighbor;
912:                     remotePoints[nL].index = remoteVertex;
913:                     ++nL;
914:                   }
915:                 }
916:               } else {
917:                 for (zv = 0; zv < nVz; ++zv) {
918:                   for (yv = 0; yv < nVy; ++yv) {
919:                     const PetscInt localVertex  = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
920:                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0     + nC; /* TODO: Correct this for neighbor sizes */

922:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
923:                       localPoints[nL]        = localVertex;
924:                       remotePoints[nL].rank  = neighbor;
925:                       remotePoints[nL].index = remoteVertex;
926:                       ++nL;
927:                     }
928:                   }
929:                 }
930: #if 0
931:                 for (xf = 0; xf < nxF; ++xf) {
932:                   /* THIS IS WRONG */
933:                   const PetscInt localFace  = 0 + nC+nV; /* right faces */
934:                   const PetscInt remoteFace = 0 + nC+nV;

936:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
937:                     localPoints[nL]        = localFace;
938:                     remotePoints[nL].rank  = neighbor;
939:                     remotePoints[nL].index = remoteFace;
940:                     ++nL;
941:                   }
942:                 }
943: #endif
944:               }
945:             }
946:           } else {
947:             if (yp < 0) { /* bottom */
948:               if (zp < 0) { /* back */
949:                 for (xv = 0; xv < nVx; ++xv) {
950:                   const PetscInt localVertex  = (      0*nVy +     0)*nVx + xv + nC; /* bottom back vertices */
951:                   const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

953:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
954:                     localPoints[nL]        = localVertex;
955:                     remotePoints[nL].rank  = neighbor;
956:                     remotePoints[nL].index = remoteVertex;
957:                     ++nL;
958:                   }
959:                 }
960:               } else if (zp > 0) { /* front */
961:                 for (xv = 0; xv < nVx; ++xv) {
962:                   const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* bottom front vertices */
963:                   const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

965:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
966:                     localPoints[nL]        = localVertex;
967:                     remotePoints[nL].rank  = neighbor;
968:                     remotePoints[nL].index = remoteVertex;
969:                     ++nL;
970:                   }
971:                 }
972:               } else {
973:                 for (zv = 0; zv < nVz; ++zv) {
974:                   for (xv = 0; xv < nVx; ++xv) {
975:                     const PetscInt localVertex  = (zv*nVy +     0)*nVx + xv + nC; /* bottom vertices */
976:                     const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

978:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
979:                       localPoints[nL]        = localVertex;
980:                       remotePoints[nL].rank  = neighbor;
981:                       remotePoints[nL].index = remoteVertex;
982:                       ++nL;
983:                     }
984:                   }
985:                 }
986: #if 0
987:                 for (yf = 0; yf < nyF; ++yf) {
988:                   /* THIS IS WRONG */
989:                   const PetscInt localFace  = 0 + nC+nV; /* bottom faces */
990:                   const PetscInt remoteFace = 0 + nC+nV;

992:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
993:                     localPoints[nL]        = localFace;
994:                     remotePoints[nL].rank  = neighbor;
995:                     remotePoints[nL].index = remoteFace;
996:                     ++nL;
997:                   }
998:                 }
999: #endif
1000:               }
1001:             } else if (yp > 0) { /* top */
1002:               if (zp < 0) { /* back */
1003:                 for (xv = 0; xv < nVx; ++xv) {
1004:                   const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
1005:                   const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

1007:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
1008:                     localPoints[nL]        = localVertex;
1009:                     remotePoints[nL].rank  = neighbor;
1010:                     remotePoints[nL].index = remoteVertex;
1011:                     ++nL;
1012:                   }
1013:                 }
1014:               } else if (zp > 0) { /* front */
1015:                 for (xv = 0; xv < nVx; ++xv) {
1016:                   const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
1017:                   const PetscInt remoteVertex = (      0*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

1019:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
1020:                     localPoints[nL]        = localVertex;
1021:                     remotePoints[nL].rank  = neighbor;
1022:                     remotePoints[nL].index = remoteVertex;
1023:                     ++nL;
1024:                   }
1025:                 }
1026:               } else {
1027:                 for (zv = 0; zv < nVz; ++zv) {
1028:                   for (xv = 0; xv < nVx; ++xv) {
1029:                     const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
1030:                     const PetscInt remoteVertex = (zv*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

1032:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
1033:                       localPoints[nL]        = localVertex;
1034:                       remotePoints[nL].rank  = neighbor;
1035:                       remotePoints[nL].index = remoteVertex;
1036:                       ++nL;
1037:                     }
1038:                   }
1039:                 }
1040: #if 0
1041:                 for (yf = 0; yf < nyF; ++yf) {
1042:                   /* THIS IS WRONG */
1043:                   const PetscInt localFace  = 0 + nC+nV; /* top faces */
1044:                   const PetscInt remoteFace = 0 + nC+nV;

1046:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1047:                     localPoints[nL]        = localFace;
1048:                     remotePoints[nL].rank  = neighbor;
1049:                     remotePoints[nL].index = remoteFace;
1050:                     ++nL;
1051:                   }
1052:                 }
1053: #endif
1054:               }
1055:             } else {
1056:               if (zp < 0) { /* back */
1057:                 for (yv = 0; yv < nVy; ++yv) {
1058:                   for (xv = 0; xv < nVx; ++xv) {
1059:                     const PetscInt localVertex  = (      0*nVy + yv)*nVx + xv + nC; /* back vertices */
1060:                     const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

1062:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
1063:                       localPoints[nL]        = localVertex;
1064:                       remotePoints[nL].rank  = neighbor;
1065:                       remotePoints[nL].index = remoteVertex;
1066:                       ++nL;
1067:                     }
1068:                   }
1069:                 }
1070: #if 0
1071:                 for (zf = 0; zf < nzF; ++zf) {
1072:                   /* THIS IS WRONG */
1073:                   const PetscInt localFace  = 0 + nC+nV; /* back faces */
1074:                   const PetscInt remoteFace = 0 + nC+nV;

1076:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1077:                     localPoints[nL]        = localFace;
1078:                     remotePoints[nL].rank  = neighbor;
1079:                     remotePoints[nL].index = remoteFace;
1080:                     ++nL;
1081:                   }
1082:                 }
1083: #endif
1084:               } else if (zp > 0) { /* front */
1085:                 for (yv = 0; yv < nVy; ++yv) {
1086:                   for (xv = 0; xv < nVx; ++xv) {
1087:                     const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
1088:                     const PetscInt remoteVertex = (      0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

1090:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
1091:                       localPoints[nL]        = localVertex;
1092:                       remotePoints[nL].rank  = neighbor;
1093:                       remotePoints[nL].index = remoteVertex;
1094:                       ++nL;
1095:                     }
1096:                   }
1097:                 }
1098: #if 0
1099:                 for (zf = 0; zf < nzF; ++zf) {
1100:                   /* THIS IS WRONG */
1101:                   const PetscInt localFace  = 0 + nC+nV; /* front faces */
1102:                   const PetscInt remoteFace = 0 + nC+nV;

1104:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1105:                     localPoints[nL]        = localFace;
1106:                     remotePoints[nL].rank  = neighbor;
1107:                     remotePoints[nL].index = remoteFace;
1108:                     ++nL;
1109:                   }
1110:                 }
1111: #endif
1112:               } else {
1113:                 /* Nothing is shared from the interior */
1114:               }
1115:             }
1116:           }
1117:         }
1118:       }
1119:     }
1120:   }
1121:   PetscBTDestroy(&isLeaf);
1122:   /* Remove duplication in leaf determination */
1123:   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);
1124:   PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);
1125:   PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1126:   DMSetPointSF(dm, sf);
1127:   PetscSFDestroy(&sf);
1128:   *s = section;
1129:   return(0);
1130: }

1132: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
1133: {
1134:   DM_DA         *da = (DM_DA *) dm->data;
1135:   Vec            coordinates;
1136:   PetscSection   section;
1137:   PetscScalar   *coords;
1138:   PetscReal      h[3];
1139:   PetscInt       dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;

1144:   DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);
1145:   if (dim > 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"The following code only works for dim <= 3");
1146:   h[0] = (xu - xl)/M;
1147:   h[1] = (yu - yl)/N;
1148:   h[2] = (zu - zl)/P;
1149:   DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
1150:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
1151:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);
1152:   PetscSectionSetNumFields(section, 1);
1153:   PetscSectionSetFieldComponents(section, 0, dim);
1154:   PetscSectionSetChart(section, vStart, vEnd);
1155:   for (v = vStart; v < vEnd; ++v) {
1156:     PetscSectionSetDof(section, v, dim);
1157:   }
1158:   PetscSectionSetUp(section);
1159:   PetscSectionGetStorageSize(section, &size);
1160:   VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
1161:   PetscObjectSetName((PetscObject)coordinates,"coordinates");
1162:   VecGetArray(coordinates, &coords);
1163:   for (k = 0; k < nVz; ++k) {
1164:     PetscInt ind[3], d, off;

1166:     ind[0] = 0;
1167:     ind[1] = 0;
1168:     ind[2] = k + da->zs;
1169:     for (j = 0; j < nVy; ++j) {
1170:       ind[1] = j + da->ys;
1171:       for (i = 0; i < nVx; ++i) {
1172:         const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;

1174:         PetscSectionGetOffset(section, vertex, &off);
1175:         ind[0] = i + da->xs;
1176:         for (d = 0; d < dim; ++d) {
1177:           coords[off+d] = h[d]*ind[d];
1178:         }
1179:       }
1180:     }
1181:   }
1182:   VecRestoreArray(coordinates, &coords);
1183:   DMSetCoordinateSection(dm, PETSC_DETERMINE, section);
1184:   DMSetCoordinatesLocal(dm, coordinates);
1185:   PetscSectionDestroy(&section);
1186:   VecDestroy(&coordinates);
1187:   return(0);
1188: }

1190: PetscErrorCode DMProjectFunctionLocal_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1191: {
1192:   PetscDS         prob;
1193:   PetscFE         fe;
1194:   PetscDualSpace  sp;
1195:   PetscQuadrature q;
1196:   PetscSection    section;
1197:   PetscScalar    *values;
1198:   PetscInt        numFields, Nc, dim, dimEmbed, spDim, totDim, numValues, cStart, cEnd, f, c, v, d;
1199:   PetscFEGeom    *geom;
1200:   PetscErrorCode  ierr;

1203:   DMGetDS(dm, &prob);
1204:   PetscDSGetTotalDimension(prob, &totDim);
1205:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1206:   PetscFEGetQuadrature(fe, &q);
1207:   DMGetSection(dm, &section);
1208:   PetscSectionGetNumFields(section, &numFields);
1209:   DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1210:   DMGetCoordinateDim(dm, &dimEmbed);
1211:   DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1212:   DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
1213:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
1214:   DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
1215:   PetscFEGeomCreate(q,1,dimEmbed,PETSC_FALSE,&geom);
1216:   for (c = cStart; c < cEnd; ++c) {
1217:     DMDAComputeCellGeometryFEM(dm, c, q, geom->v, geom->J, NULL, geom->detJ);
1218:     for (f = 0, v = 0; f < numFields; ++f) {
1219:       void * const ctx = ctxs ? ctxs[f] : NULL;

1221:       PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1222:       PetscFEGetNumComponents(fe, &Nc);
1223:       PetscFEGetDualSpace(fe, &sp);
1224:       PetscDualSpaceGetDimension(sp, &spDim);
1225:       for (d = 0; d < spDim; ++d) {
1226:         PetscDualSpaceApply(sp, d, time, geom, Nc, funcs[f], ctx, &values[v]);
1227:       }
1228:     }
1229:     DMDAVecSetClosure(dm, section, localX, c, values, mode);
1230:   }
1231:   PetscFEGeomDestroy(&geom);
1232:   DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
1233:   return(0);
1234: }

1236: PetscErrorCode DMComputeL2Diff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1237: {
1238:   const PetscInt  debug = 0;
1239:   PetscDS    prob;
1240:   PetscFE         fe;
1241:   PetscQuadrature quad;
1242:   PetscSection    section;
1243:   Vec             localX;
1244:   PetscScalar    *funcVal;
1245:   PetscReal      *coords, *v0, *J, *invJ, detJ;
1246:   PetscReal       localDiff = 0.0;
1247:   PetscInt        dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1248:   PetscErrorCode  ierr;

1251:   DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1252:   DMGetDS(dm, &prob);
1253:   PetscDSGetTotalComponents(prob, &Nc);
1254:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1255:   PetscFEGetQuadrature(fe, &quad);
1256:   DMGetSection(dm, &section);
1257:   PetscSectionGetNumFields(section, &numFields);
1258:   DMGetLocalVector(dm, &localX);
1259:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1260:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1261:   /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1262:   PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1263:   DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1264:   for (c = cStart; c < cEnd; ++c) {
1265:     PetscScalar *x = NULL;
1266:     PetscReal    elemDiff = 0.0;

1268:     DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);
1269:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1270:     DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);

1272:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1273:       void * const ctx = ctxs ? ctxs[field] : NULL;
1274:       const PetscReal *quadPoints, *quadWeights;
1275:       PetscReal       *basis;
1276:       PetscInt         qNc, Nq, Nb, Nc, q, d, e, fc, f;

1278:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1279:       PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1280:       PetscFEGetDimension(fe, &Nb);
1281:       PetscFEGetNumComponents(fe, &Nc);
1282:       if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Quadrature components %D != %D field components\n", qNc, Nc);
1283:       PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
1284:       if (debug) {
1285:         char title[1024];
1286:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1287:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1288:       }
1289:       for (q = 0; q < Nq; ++q) {
1290:         for (d = 0; d < dim; d++) {
1291:           coords[d] = v0[d];
1292:           for (e = 0; e < dim; e++) {
1293:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1294:           }
1295:         }
1296:         (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);
1297:         for (fc = 0; fc < Nc; ++fc) {
1298:           PetscScalar interpolant = 0.0;

1300:           for (f = 0; f < Nb; ++f) {
1301:             interpolant += x[fieldOffset+f]*basis[(q*Nb+f)*Nc+fc];
1302:           }
1303:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+fc]*detJ);}
1304:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ;
1305:         }
1306:       }
1307:       comp        += Nc;
1308:       fieldOffset += Nb;
1309:     }
1310:     DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1311:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
1312:     localDiff += elemDiff;
1313:   }
1314:   PetscFree5(funcVal,coords,v0,J,invJ);
1315:   DMRestoreLocalVector(dm, &localX);
1316:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1317:   *diff = PetscSqrtReal(*diff);
1318:   return(0);
1319: }

1321: PetscErrorCode DMComputeL2GradientDiff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
1322: {
1323:   const PetscInt  debug = 0;
1324:   PetscDS    prob;
1325:   PetscFE         fe;
1326:   PetscQuadrature quad;
1327:   PetscSection    section;
1328:   Vec             localX;
1329:   PetscScalar    *funcVal, *interpolantVec;
1330:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
1331:   PetscReal       localDiff = 0.0;
1332:   PetscInt        dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1333:   PetscErrorCode  ierr;

1336:   DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1337:   DMGetDS(dm, &prob);
1338:   PetscDSGetTotalComponents(prob, &Nc);
1339:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1340:   PetscFEGetQuadrature(fe, &quad);
1341:   DMGetSection(dm, &section);
1342:   PetscSectionGetNumFields(section, &numFields);
1343:   DMGetLocalVector(dm, &localX);
1344:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1345:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1346:   /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1347:   PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
1348:   DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1349:   for (c = cStart; c < cEnd; ++c) {
1350:     PetscScalar *x = NULL;
1351:     PetscReal    elemDiff = 0.0;

1353:     DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);
1354:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1355:     DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);

1357:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1358:       void * const ctx = ctxs ? ctxs[field] : NULL;
1359:       const PetscReal *quadPoints, *quadWeights;
1360:       PetscReal       *basisDer;
1361:       PetscInt         qNc, Nq, Nb, Nc, q, d, e, fc, f, g;

1363:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1364:       PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1365:       PetscFEGetDimension(fe, &Nb);
1366:       PetscFEGetNumComponents(fe, &Nc);
1367:       if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Quadrature components %D != %D field components\n", qNc, Nc);
1368:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1369:       if (debug) {
1370:         char title[1024];
1371:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1372:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1373:       }
1374:       for (q = 0; q < Nq; ++q) {
1375:         for (d = 0; d < dim; d++) {
1376:           coords[d] = v0[d];
1377:           for (e = 0; e < dim; e++) {
1378:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1379:           }
1380:         }
1381:         (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);
1382:         for (fc = 0; fc < Nc; ++fc) {
1383:           PetscScalar interpolant = 0.0;

1385:           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
1386:           for (f = 0; f < Nb; ++f) {
1387:             for (d = 0; d < dim; ++d) {
1388:               realSpaceDer[d] = 0.0;
1389:               for (g = 0; g < dim; ++g) {
1390:                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb+f)*dim+g];
1391:               }
1392:               interpolantVec[d] += x[fieldOffset+f]*realSpaceDer[d];
1393:             }
1394:           }
1395:           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
1396:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ);}
1397:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ;
1398:         }
1399:       }
1400:       comp        += Nc;
1401:       fieldOffset += Nb*Nc;
1402:     }
1403:     DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1404:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
1405:     localDiff += elemDiff;
1406:   }
1407:   PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
1408:   DMRestoreLocalVector(dm, &localX);
1409:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1410:   *diff = PetscSqrtReal(*diff);
1411:   return(0);
1412: }

1414: /* ------------------------------------------------------------------- */

1416: /*@C
1417:      DMDAGetArray - Gets a work array for a DMDA

1419:     Input Parameter:
1420: +    da - information about my local patch
1421: -    ghosted - do you want arrays for the ghosted or nonghosted patch

1423:     Output Parameters:
1424: .    vptr - array data structured

1426:     Note:  The vector values are NOT initialized and may have garbage in them, so you may need
1427:            to zero them.

1429:   Level: advanced

1431: .seealso: DMDARestoreArray()

1433: @*/
1434: PetscErrorCode  DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
1435: {
1437:   PetscInt       j,i,xs,ys,xm,ym,zs,zm;
1438:   char           *iarray_start;
1439:   void           **iptr = (void**)vptr;
1440:   DM_DA          *dd    = (DM_DA*)da->data;

1444:   if (ghosted) {
1445:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1446:       if (dd->arrayghostedin[i]) {
1447:         *iptr                 = dd->arrayghostedin[i];
1448:         iarray_start          = (char*)dd->startghostedin[i];
1449:         dd->arrayghostedin[i] = NULL;
1450:         dd->startghostedin[i] = NULL;

1452:         goto done;
1453:       }
1454:     }
1455:     xs = dd->Xs;
1456:     ys = dd->Ys;
1457:     zs = dd->Zs;
1458:     xm = dd->Xe-dd->Xs;
1459:     ym = dd->Ye-dd->Ys;
1460:     zm = dd->Ze-dd->Zs;
1461:   } else {
1462:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1463:       if (dd->arrayin[i]) {
1464:         *iptr          = dd->arrayin[i];
1465:         iarray_start   = (char*)dd->startin[i];
1466:         dd->arrayin[i] = NULL;
1467:         dd->startin[i] = NULL;

1469:         goto done;
1470:       }
1471:     }
1472:     xs = dd->xs;
1473:     ys = dd->ys;
1474:     zs = dd->zs;
1475:     xm = dd->xe-dd->xs;
1476:     ym = dd->ye-dd->ys;
1477:     zm = dd->ze-dd->zs;
1478:   }

1480:   switch (da->dim) {
1481:   case 1: {
1482:     void *ptr;

1484:     PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);

1486:     ptr   = (void*)(iarray_start - xs*sizeof(PetscScalar));
1487:     *iptr = (void*)ptr;
1488:     break;
1489:   }
1490:   case 2: {
1491:     void **ptr;

1493:     PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);

1495:     ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
1496:     for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
1497:     *iptr = (void*)ptr;
1498:     break;
1499:   }
1500:   case 3: {
1501:     void ***ptr,**bptr;

1503:     PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);

1505:     ptr  = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1506:     bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1507:     for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1508:     for (i=zs; i<zs+zm; i++) {
1509:       for (j=ys; j<ys+ym; j++) {
1510:         ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1511:       }
1512:     }
1513:     *iptr = (void*)ptr;
1514:     break;
1515:   }
1516:   default:
1517:     SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
1518:   }

1520: done:
1521:   /* add arrays to the checked out list */
1522:   if (ghosted) {
1523:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1524:       if (!dd->arrayghostedout[i]) {
1525:         dd->arrayghostedout[i] = *iptr;
1526:         dd->startghostedout[i] = iarray_start;
1527:         break;
1528:       }
1529:     }
1530:   } else {
1531:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1532:       if (!dd->arrayout[i]) {
1533:         dd->arrayout[i] = *iptr;
1534:         dd->startout[i] = iarray_start;
1535:         break;
1536:       }
1537:     }
1538:   }
1539:   return(0);
1540: }

1542: /*@C
1543:      DMDARestoreArray - Restores an array of derivative types for a DMDA

1545:     Input Parameter:
1546: +    da - information about my local patch
1547: .    ghosted - do you want arrays for the ghosted or nonghosted patch
1548: -    vptr - array data structured to be passed to ad_FormFunctionLocal()

1550:      Level: advanced

1552: .seealso: DMDAGetArray()

1554: @*/
1555: PetscErrorCode  DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1556: {
1557:   PetscInt i;
1558:   void     **iptr = (void**)vptr,*iarray_start = 0;
1559:   DM_DA    *dd    = (DM_DA*)da->data;

1563:   if (ghosted) {
1564:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1565:       if (dd->arrayghostedout[i] == *iptr) {
1566:         iarray_start           = dd->startghostedout[i];
1567:         dd->arrayghostedout[i] = NULL;
1568:         dd->startghostedout[i] = NULL;
1569:         break;
1570:       }
1571:     }
1572:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1573:       if (!dd->arrayghostedin[i]) {
1574:         dd->arrayghostedin[i] = *iptr;
1575:         dd->startghostedin[i] = iarray_start;
1576:         break;
1577:       }
1578:     }
1579:   } else {
1580:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1581:       if (dd->arrayout[i] == *iptr) {
1582:         iarray_start    = dd->startout[i];
1583:         dd->arrayout[i] = NULL;
1584:         dd->startout[i] = NULL;
1585:         break;
1586:       }
1587:     }
1588:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1589:       if (!dd->arrayin[i]) {
1590:         dd->arrayin[i] = *iptr;
1591:         dd->startin[i] = iarray_start;
1592:         break;
1593:       }
1594:     }
1595:   }
1596:   return(0);
1597: }