Actual source code: dalocal.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: /*
  3:   Code for manipulating distributed regular arrays in parallel.
  4: */

  6: #include <petsc-private/dmdaimpl.h>    /*I   "petscdmda.h"   I*/
  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 */
 20: static PetscErrorCode  VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
 21: {
 23:   PetscInt       n,m;
 24:   Vec            vec = (Vec)obj;
 25:   PetscScalar    *array;
 26:   mxArray        *mat;
 27:   DM             da;

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

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

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


 52: PetscErrorCode  DMCreateLocalVector_DA(DM da,Vec *g)
 53: {
 55:   DM_DA          *dd = (DM_DA*)da->data;

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

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

 82:   Input Parameter:
 83: . dm - The DM object

 85:   Output Parameters:
 86: + numCellsX - The number of local cells in the x-direction
 87: . numCellsY - The number of local cells in the y-direction
 88: . numCellsZ - The number of local cells in the z-direction
 89: - numCells - The number of local cells

 91:   Level: developer

 93: .seealso: DMDAGetCellPoint()
 94: @*/
 95: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
 96: {
 97:   DM_DA         *da  = (DM_DA*) dm->data;
 98:   const PetscInt dim = da->dim;
 99:   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
100:   const PetscInt nC  = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);

104:   if (numCellsX) {
106:     *numCellsX = mx;
107:   }
108:   if (numCellsY) {
110:     *numCellsY = my;
111:   }
112:   if (numCellsZ) {
114:     *numCellsZ = mz;
115:   }
116:   if (numCells) {
118:     *numCells = nC;
119:   }
120:   return(0);
121: }

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

128:   Input Parameters:
129: + dm - The DM object
130: - i,j,k - The global indices for the cell

132:   Output Parameters:
133: . point - The local DM point

135:   Level: developer

137: .seealso: DMDAGetNumCells()
138: @*/
139: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
140: {
141:   DM_DA         *da  = (DM_DA*) dm->data;
142:   const PetscInt dim = da->dim;
143:   DMDALocalInfo  info;

149:   DMDAGetLocalInfo(dm, &info);
150:   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);}
151:   if (dim > 1) {if ((i < info.gys) || (i >= info.gys+info.gym)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %d not in [%d, %d)", i, info.gys, info.gys+info.gym);}
152:   if (dim > 2) {if ((i < info.gzs) || (i >= info.gzs+info.gzm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %d not in [%d, %d)", i, info.gzs, info.gzs+info.gzm);}
153:   *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0);
154:   return(0);
155: }

159: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
160: {
161:   DM_DA          *da = (DM_DA*) dm->data;
162:   const PetscInt dim = da->dim;
163:   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
164:   const PetscInt nVx = mx+1;
165:   const PetscInt nVy = dim > 1 ? (my+1) : 1;
166:   const PetscInt nVz = dim > 2 ? (mz+1) : 1;
167:   const PetscInt nV  = nVx*nVy*nVz;

170:   if (numVerticesX) {
172:     *numVerticesX = nVx;
173:   }
174:   if (numVerticesY) {
176:     *numVerticesY = nVy;
177:   }
178:   if (numVerticesZ) {
180:     *numVerticesZ = nVz;
181:   }
182:   if (numVertices) {
184:     *numVertices = nV;
185:   }
186:   return(0);
187: }

191: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
192: {
193:   DM_DA          *da = (DM_DA*) dm->data;
194:   const PetscInt dim = da->dim;
195:   const PetscInt mx  = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
196:   const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
197:   const PetscInt nXF = (mx+1)*nxF;
198:   const PetscInt nyF = mx*(dim > 2 ? mz : 1);
199:   const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0;
200:   const PetscInt nzF = mx*(dim > 1 ? my : 0);
201:   const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0;

204:   if (numXFacesX) {
206:     *numXFacesX = nxF;
207:   }
208:   if (numXFaces) {
210:     *numXFaces = nXF;
211:   }
212:   if (numYFacesY) {
214:     *numYFacesY = nyF;
215:   }
216:   if (numYFaces) {
218:     *numYFaces = nYF;
219:   }
220:   if (numZFacesZ) {
222:     *numZFacesZ = nzF;
223:   }
224:   if (numZFaces) {
226:     *numZFaces = nZF;
227:   }
228:   return(0);
229: }

233: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
234: {
235:   DM_DA          *da = (DM_DA*) dm->data;
236:   const PetscInt dim = da->dim;
237:   PetscInt       nC, nV, nXF, nYF, nZF;

243:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
244:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
245:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
246:   if (height == 0) {
247:     /* Cells */
248:     if (pStart) *pStart = 0;
249:     if (pEnd)   *pEnd   = nC;
250:   } else if (height == 1) {
251:     /* Faces */
252:     if (pStart) *pStart = nC+nV;
253:     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
254:   } else if (height == dim) {
255:     /* Vertices */
256:     if (pStart) *pStart = nC;
257:     if (pEnd)   *pEnd   = nC+nV;
258:   } else if (height < 0) {
259:     /* All points */
260:     if (pStart) *pStart = 0;
261:     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
262:   } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
263:   return(0);
264: }

268: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
269: {
270:   DM_DA         *da  = (DM_DA*) dm->data;
271:   const PetscInt dim = da->dim;
272:   PetscInt       nC, nV, nXF, nYF, nZF;

278:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
279:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
280:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
281:   if (depth == dim) {
282:     /* Cells */
283:     if (pStart) *pStart = 0;
284:     if (pEnd)   *pEnd   = nC;
285:   } else if (depth == dim-1) {
286:     /* Faces */
287:     if (pStart) *pStart = nC+nV;
288:     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
289:   } else if (depth == 0) {
290:     /* Vertices */
291:     if (pStart) *pStart = nC;
292:     if (pEnd)   *pEnd   = nC+nV;
293:   } else if (depth < 0) {
294:     /* All points */
295:     if (pStart) *pStart = 0;
296:     if (pEnd)   *pEnd   = nC+nV+nXF+nYF+nZF;
297:   } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
298:   return(0);
299: }

303: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
304: {
305:   DM_DA         *da  = (DM_DA*) dm->data;
306:   const PetscInt dim = da->dim;
307:   PetscInt       nC, nV, nXF, nYF, nZF;

311:   *coneSize = 0;
312:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
313:   DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
314:   DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
315:   switch (dim) {
316:   case 2:
317:     if (p >= 0) {
318:       if (p < nC) {
319:         *coneSize = 4;
320:       } else if (p < nC+nV) {
321:         *coneSize = 0;
322:       } else if (p < nC+nV+nXF+nYF+nZF) {
323:         *coneSize = 2;
324:       } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
325:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
326:     break;
327:   case 3:
328:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
329:     break;
330:   }
331:   return(0);
332: }

336: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
337: {
338:   DM_DA         *da  = (DM_DA*) dm->data;
339:   const PetscInt dim = da->dim;
340:   PetscInt       nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;

344:   if (!cone) {DMGetWorkArray(dm, 6, PETSC_INT, cone);}
345:   DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);
346:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
347:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
348:   switch (dim) {
349:   case 2:
350:     if (p >= 0) {
351:       if (p < nC) {
352:         const PetscInt cy = p / nCx;
353:         const PetscInt cx = p % nCx;

355:         (*cone)[0] = cy*nxF + cx + nC+nV;
356:         (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
357:         (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
358:         (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
359:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
360:       } else if (p < nC+nV) {
361:       } else if (p < nC+nV+nXF) {
362:         const PetscInt fy = (p - nC+nV) / nxF;
363:         const PetscInt fx = (p - nC+nV) % nxF;

365:         (*cone)[0] = fy*nVx + fx + nC;
366:         (*cone)[1] = fy*nVx + fx + 1 + nC;
367:       } else if (p < nC+nV+nXF+nYF) {
368:         const PetscInt fx = (p - nC+nV+nXF) / nyF;
369:         const PetscInt fy = (p - nC+nV+nXF) % nyF;

371:         (*cone)[0] = fy*nVx + fx + nC;
372:         (*cone)[1] = fy*nVx + fx + nVx + nC;
373:       } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
374:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
375:     break;
376:   case 3:
377:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
378:     break;
379:   }
380:   return(0);
381: }

385: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
386: {

390:   DMGetWorkArray(dm, 6, PETSC_INT, cone);
391:   return(0);
392: }

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

400:   Input Parameters:
401: + dm- The DMDA
402: . numFields - The number of fields
403: . numComp - The number of components in each field
404: . numDof - The number of dofs per dimension for each field
405: . numFaceDof - The number of dofs per face for each field and direction, or NULL

407:   Level: developer

409:   Note:
410:   The default DMDA numbering is as follows:

412:     - Cells:    [0,             nC)
413:     - Vertices: [nC,            nC+nV)
414:     - X-Faces:  [nC+nV,         nC+nV+nXF)         normal is +- x-dir
415:     - Y-Faces:  [nC+nV+nXF,     nC+nV+nXF+nYF)     normal is +- y-dir
416:     - Z-Faces:  [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir

418:   We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment.
419: @*/
420: PetscErrorCode DMDACreateSection(DM dm, const PetscInt numComp[], const PetscInt numDof[], const PetscInt numFaceDof[], PetscSection *s)
421: {
422:   DM_DA            *da  = (DM_DA*) dm->data;
423:   PetscSection      section;
424:   const PetscInt    dim = da->dim;
425:   PetscInt          numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0};
426:   PetscBT           isLeaf;
427:   PetscSF           sf;
428:   PetscMPIInt       rank;
429:   const PetscMPIInt *neighbors;
430:   PetscInt          *localPoints;
431:   PetscSFNode       *remotePoints;
432:   PetscInt          nleaves = 0,  nleavesCheck = 0, nL = 0;
433:   PetscInt          nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF;
434:   PetscInt          pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd;
435:   PetscInt          f, v, c, xf, yf, zf, xn, yn, zn;
436:   PetscErrorCode    ierr;

441:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
442:   DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
443:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
444:   DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
445:   DMDAGetHeightStratum(dm, -1,  &pStart, &pEnd);
446:   DMDAGetHeightStratum(dm, 0,   &cStart, &cEnd);
447:   DMDAGetHeightStratum(dm, 1,   &fStart, &fEnd);
448:   DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
449:   xfStart = vEnd;  xfEnd = xfStart+nXF;
450:   yfStart = xfEnd; yfEnd = yfStart+nYF;
451:   zfStart = yfEnd; zfEnd = zfStart+nZF;
452:   if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd);
453:   /* Create local section */
454:   DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);
455:   for (f = 0; f < numFields; ++f) {
456:     numVertexTotDof  += numDof[f*(dim+1)+0];
457:     numCellTotDof    += numDof[f*(dim+1)+dim];
458:     numFaceTotDof[0] += dim > 0 ? (numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]) : 0;
459:     numFaceTotDof[1] += dim > 1 ? (numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]) : 0;
460:     numFaceTotDof[2] += dim > 2 ? (numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]) : 0;
461:   }
462:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
463:   if (numFields > 0) {
464:     PetscSectionSetNumFields(section, numFields);
465:     for (f = 0; f < numFields; ++f) {
466:       const char *name;

468:       DMDAGetFieldName(dm, f, &name);
469:       PetscSectionSetFieldName(section, f, name ? name : "Field");
470:       if (numComp) {
471:         PetscSectionSetFieldComponents(section, f, numComp[f]);
472:       }
473:     }
474:   }
475:   PetscSectionSetChart(section, pStart, pEnd);
476:   for (v = vStart; v < vEnd; ++v) {
477:     for (f = 0; f < numFields; ++f) {
478:       PetscSectionSetFieldDof(section, v, f, numDof[f*(dim+1)+0]);
479:     }
480:     PetscSectionSetDof(section, v, numVertexTotDof);
481:   }
482:   for (xf = xfStart; xf < xfEnd; ++xf) {
483:     for (f = 0; f < numFields; ++f) {
484:       PetscSectionSetFieldDof(section, xf, f, numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]);
485:     }
486:     PetscSectionSetDof(section, xf, numFaceTotDof[0]);
487:   }
488:   for (yf = yfStart; yf < yfEnd; ++yf) {
489:     for (f = 0; f < numFields; ++f) {
490:       PetscSectionSetFieldDof(section, yf, f, numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]);
491:     }
492:     PetscSectionSetDof(section, yf, numFaceTotDof[1]);
493:   }
494:   for (zf = zfStart; zf < zfEnd; ++zf) {
495:     for (f = 0; f < numFields; ++f) {
496:       PetscSectionSetFieldDof(section, zf, f, numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]);
497:     }
498:     PetscSectionSetDof(section, zf, numFaceTotDof[2]);
499:   }
500:   for (c = cStart; c < cEnd; ++c) {
501:     for (f = 0; f < numFields; ++f) {
502:       PetscSectionSetFieldDof(section, c, f, numDof[f*(dim+1)+dim]);
503:     }
504:     PetscSectionSetDof(section, c, numCellTotDof);
505:   }
506:   PetscSectionSetUp(section);
507:   /* Create mesh point SF */
508:   PetscBTCreate(pEnd-pStart, &isLeaf);
509:   DMDAGetNeighbors(dm, &neighbors);
510:   for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
511:     for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
512:       for (xn = 0; xn < 3; ++xn) {
513:         const PetscInt xp       = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
514:         const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
515:         PetscInt       xv, yv, zv;

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

726:         if (neighbor >= 0 && neighbor < rank) {
727:           if (xp < 0) { /* left */
728:             if (yp < 0) { /* bottom */
729:               if (zp < 0) { /* back */
730:                 const PetscInt localVertex  = (      0*nVy +     0)*nVx +     0 + nC; /* left bottom back vertex */
731:                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

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

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

754:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
755:                     localPoints[nL]        = localVertex;
756:                     remotePoints[nL].rank  = neighbor;
757:                     remotePoints[nL].index = remoteVertex;
758:                     ++nL;
759:                   }
760:                 }
761:               }
762:             } else if (yp > 0) { /* top */
763:               if (zp < 0) { /* back */
764:                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx +     0 + nC; /* left top back vertex */
765:                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

767:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
768:                   localPoints[nL]        = localVertex;
769:                   remotePoints[nL].rank  = neighbor;
770:                   remotePoints[nL].index = remoteVertex;
771:                   ++nL;
772:                 }
773:               } else if (zp > 0) { /* front */
774:                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* left top front vertex */
775:                 const PetscInt remoteVertex = (      0*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

777:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
778:                   localPoints[nL]        = localVertex;
779:                   remotePoints[nL].rank  = neighbor;
780:                   remotePoints[nL].index = remoteVertex;
781:                   ++nL;
782:                 }
783:               } else {
784:                 for (zv = 0; zv < nVz; ++zv) {
785:                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx +     0 + nC; /* left top vertices */
786:                   const PetscInt remoteVertex = (zv*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

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

802:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
803:                     localPoints[nL]        = localVertex;
804:                     remotePoints[nL].rank  = neighbor;
805:                     remotePoints[nL].index = remoteVertex;
806:                     ++nL;
807:                   }
808:                 }
809:               } else if (zp > 0) { /* front */
810:                 for (yv = 0; yv < nVy; ++yv) {
811:                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx +     0 + nC; /* left front vertices */
812:                   const PetscInt remoteVertex = (      0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

814:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
815:                     localPoints[nL]        = localVertex;
816:                     remotePoints[nL].rank  = neighbor;
817:                     remotePoints[nL].index = remoteVertex;
818:                     ++nL;
819:                   }
820:                 }
821:               } else {
822:                 for (zv = 0; zv < nVz; ++zv) {
823:                   for (yv = 0; yv < nVy; ++yv) {
824:                     const PetscInt localVertex  = (zv*nVy + yv)*nVx +     0 + nC; /* left vertices */
825:                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + 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:                   }
834:                 }
835: #if 0
836:                 for (xf = 0; xf < nxF; ++xf) {
837:                   /* THIS IS WRONG */
838:                   const PetscInt localFace  = 0 + nC+nV; /* left faces */
839:                   const PetscInt remoteFace = 0 + nC+nV;

841:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
842:                     localPoints[nL]        = localFace;
843:                     remotePoints[nL].rank  = neighbor;
844:                     remotePoints[nL].index = remoteFace;
845:                   }
846:                 }
847: #endif
848:               }
849:             }
850:           } else if (xp > 0) { /* right */
851:             if (yp < 0) { /* bottom */
852:               if (zp < 0) { /* back */
853:                 const PetscInt localVertex  = (      0*nVy +     0)*nVx + nVx-1 + nC; /* right bottom back vertex */
854:                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

856:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
857:                   localPoints[nL]        = localVertex;
858:                   remotePoints[nL].rank  = neighbor;
859:                   remotePoints[nL].index = remoteVertex;
860:                   ++nL;
861:                 }
862:               } else if (zp > 0) { /* front */
863:                 const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* right bottom front vertex */
864:                 const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

866:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
867:                   localPoints[nL]        = localVertex;
868:                   remotePoints[nL].rank  = neighbor;
869:                   remotePoints[nL].index = remoteVertex;
870:                   ++nL;
871:                 }
872:               } else {
873:                 nleavesCheck += nVz;
874:                 for (zv = 0; zv < nVz; ++zv) {
875:                   const PetscInt localVertex  = (zv*nVy +     0)*nVx + nVx-1 + nC; /* right bottom vertices */
876:                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

878:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
879:                     localPoints[nL]        = localVertex;
880:                     remotePoints[nL].rank  = neighbor;
881:                     remotePoints[nL].index = remoteVertex;
882:                     ++nL;
883:                   }
884:                 }
885:               }
886:             } else if (yp > 0) { /* top */
887:               if (zp < 0) { /* back */
888:                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
889:                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

891:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
892:                   localPoints[nL]        = localVertex;
893:                   remotePoints[nL].rank  = neighbor;
894:                   remotePoints[nL].index = remoteVertex;
895:                   ++nL;
896:                 }
897:               } else if (zp > 0) { /* front */
898:                 const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
899:                 const PetscInt remoteVertex = (      0*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

901:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
902:                   localPoints[nL]        = localVertex;
903:                   remotePoints[nL].rank  = neighbor;
904:                   remotePoints[nL].index = remoteVertex;
905:                   ++nL;
906:                 }
907:               } else {
908:                 for (zv = 0; zv < nVz; ++zv) {
909:                   const PetscInt localVertex  = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
910:                   const PetscInt remoteVertex = (zv*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

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

926:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
927:                     localPoints[nL]        = localVertex;
928:                     remotePoints[nL].rank  = neighbor;
929:                     remotePoints[nL].index = remoteVertex;
930:                     ++nL;
931:                   }
932:                 }
933:               } else if (zp > 0) { /* front */
934:                 for (yv = 0; yv < nVy; ++yv) {
935:                   const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
936:                   const PetscInt remoteVertex = (      0*nVy + yv)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

938:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
939:                     localPoints[nL]        = localVertex;
940:                     remotePoints[nL].rank  = neighbor;
941:                     remotePoints[nL].index = remoteVertex;
942:                     ++nL;
943:                   }
944:                 }
945:               } else {
946:                 for (zv = 0; zv < nVz; ++zv) {
947:                   for (yv = 0; yv < nVy; ++yv) {
948:                     const PetscInt localVertex  = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
949:                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0     + nC; /* TODO: Correct this for neighbor sizes */

951:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
952:                       localPoints[nL]        = localVertex;
953:                       remotePoints[nL].rank  = neighbor;
954:                       remotePoints[nL].index = remoteVertex;
955:                       ++nL;
956:                     }
957:                   }
958:                 }
959: #if 0
960:                 for (xf = 0; xf < nxF; ++xf) {
961:                   /* THIS IS WRONG */
962:                   const PetscInt localFace  = 0 + nC+nV; /* right faces */
963:                   const PetscInt remoteFace = 0 + nC+nV;

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

982:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
983:                     localPoints[nL]        = localVertex;
984:                     remotePoints[nL].rank  = neighbor;
985:                     remotePoints[nL].index = remoteVertex;
986:                     ++nL;
987:                   }
988:                 }
989:               } else if (zp > 0) { /* front */
990:                 for (xv = 0; xv < nVx; ++xv) {
991:                   const PetscInt localVertex  = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* bottom front vertices */
992:                   const PetscInt remoteVertex = (      0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

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 {
1002:                 for (zv = 0; zv < nVz; ++zv) {
1003:                   for (xv = 0; xv < nVx; ++xv) {
1004:                     const PetscInt localVertex  = (zv*nVy +     0)*nVx + xv + nC; /* bottom vertices */
1005:                     const PetscInt remoteVertex = (zv*nVy + nVy-1)*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:                 }
1015: #if 0
1016:                 for (yf = 0; yf < nyF; ++yf) {
1017:                   /* THIS IS WRONG */
1018:                   const PetscInt localFace  = 0 + nC+nV; /* bottom faces */
1019:                   const PetscInt remoteFace = 0 + nC+nV;

1021:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1022:                     localPoints[nL]        = localFace;
1023:                     remotePoints[nL].rank  = neighbor;
1024:                     remotePoints[nL].index = remoteFace;
1025:                     ++nL;
1026:                   }
1027:                 }
1028: #endif
1029:               }
1030:             } else if (yp > 0) { /* top */
1031:               if (zp < 0) { /* back */
1032:                 for (xv = 0; xv < nVx; ++xv) {
1033:                   const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
1034:                   const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

1036:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
1037:                     localPoints[nL]        = localVertex;
1038:                     remotePoints[nL].rank  = neighbor;
1039:                     remotePoints[nL].index = remoteVertex;
1040:                     ++nL;
1041:                   }
1042:                 }
1043:               } else if (zp > 0) { /* front */
1044:                 for (xv = 0; xv < nVx; ++xv) {
1045:                   const PetscInt localVertex  = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
1046:                   const PetscInt remoteVertex = (      0*nVy +     0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

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

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

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

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

1105:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1106:                     localPoints[nL]        = localFace;
1107:                     remotePoints[nL].rank  = neighbor;
1108:                     remotePoints[nL].index = remoteFace;
1109:                     ++nL;
1110:                   }
1111:                 }
1112: #endif
1113:               } else if (zp > 0) { /* front */
1114:                 for (yv = 0; yv < nVy; ++yv) {
1115:                   for (xv = 0; xv < nVx; ++xv) {
1116:                     const PetscInt localVertex  = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
1117:                     const PetscInt remoteVertex = (      0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

1119:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
1120:                       localPoints[nL]        = localVertex;
1121:                       remotePoints[nL].rank  = neighbor;
1122:                       remotePoints[nL].index = remoteVertex;
1123:                       ++nL;
1124:                     }
1125:                   }
1126:                 }
1127: #if 0
1128:                 for (zf = 0; zf < nzF; ++zf) {
1129:                   /* THIS IS WRONG */
1130:                   const PetscInt localFace  = 0 + nC+nV; /* front faces */
1131:                   const PetscInt remoteFace = 0 + nC+nV;

1133:                   if (!PetscBTLookupSet(isLeaf, localFace)) {
1134:                     localPoints[nL]        = localFace;
1135:                     remotePoints[nL].rank  = neighbor;
1136:                     remotePoints[nL].index = remoteFace;
1137:                     ++nL;
1138:                   }
1139:                 }
1140: #endif
1141:               } else {
1142:                 /* Nothing is shared from the interior */
1143:               }
1144:             }
1145:           }
1146:         }
1147:       }
1148:     }
1149:   }
1150:   PetscBTDestroy(&isLeaf);
1151:   /* Remove duplication in leaf determination */
1152:   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);
1153:   PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);
1154:   PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1155:   DMSetPointSF(dm, sf);
1156:   PetscSFDestroy(&sf);
1157:   *s = section;
1158:   return(0);
1159: }

1163: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
1164: {
1165:   DM_DA         *da = (DM_DA *) dm->data;
1166:   Vec            coordinates;
1167:   PetscSection   section;
1168:   PetscScalar   *coords;
1169:   PetscReal      h[3];
1170:   PetscInt       dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;

1175:   DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);
1176:   h[0] = (xu - xl)/M;
1177:   h[1] = (yu - yl)/N;
1178:   h[2] = (zu - zl)/P;
1179:   DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
1180:   DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
1181:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);
1182:   PetscSectionSetNumFields(section, 1);
1183:   PetscSectionSetFieldComponents(section, 0, dim);
1184:   PetscSectionSetChart(section, vStart, vEnd);
1185:   for (v = vStart; v < vEnd; ++v) {
1186:     PetscSectionSetDof(section, v, dim);
1187:   }
1188:   PetscSectionSetUp(section);
1189:   PetscSectionGetStorageSize(section, &size);
1190:   VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
1191:   VecGetArray(coordinates, &coords);
1192:   for (k = 0; k < nVz; ++k) {
1193:     PetscInt ind[3], d, off;

1195:     ind[0] = 0;
1196:     ind[1] = 0;
1197:     ind[2] = k + da->zs;
1198:     for (j = 0; j < nVy; ++j) {
1199:       ind[1] = j + da->ys;
1200:       for (i = 0; i < nVx; ++i) {
1201:         const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;

1203:         PetscSectionGetOffset(section, vertex, &off);
1204:         ind[0] = i + da->xs;
1205:         for (d = 0; d < dim; ++d) {
1206:           coords[off+d] = h[d]*ind[d];
1207:         }
1208:       }
1209:     }
1210:   }
1211:   VecRestoreArray(coordinates, &coords);
1212:   DMSetCoordinateSection(dm, section);
1213:   DMSetCoordinatesLocal(dm, coordinates);
1214:   PetscSectionDestroy(&section);
1215:   VecDestroy(&coordinates);
1216:   return(0);
1217: }

1221: PetscErrorCode DMDAProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1222: {
1223:   PetscDS    prob;
1224:   PetscFE         fe;
1225:   PetscDualSpace  sp;
1226:   PetscQuadrature q;
1227:   PetscSection    section;
1228:   PetscScalar    *values;
1229:   PetscReal      *v0, *J, *detJ;
1230:   PetscInt        numFields, numComp, numPoints, dim, spDim, totDim, numValues, cStart, cEnd, f, c, v, d;
1231:   PetscErrorCode  ierr;

1234:   DMGetDS(dm, &prob);
1235:   PetscDSGetTotalDimension(prob, &totDim);
1236:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1237:   PetscFEGetQuadrature(fe, &q);
1238:   DMGetDefaultSection(dm, &section);
1239:   PetscSectionGetNumFields(section, &numFields);
1240:   DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1241:   DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1242:   DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
1243:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
1244:   DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
1245:   PetscQuadratureGetData(q, NULL, &numPoints, NULL, NULL);
1246:   PetscMalloc3(dim*numPoints,&v0,dim*dim*numPoints,&J,numPoints,&detJ);
1247:   for (c = cStart; c < cEnd; ++c) {
1248:     PetscCellGeometry geom;

1250:     DMDAComputeCellGeometry(dm, c, q, v0, J, NULL, detJ);
1251:     geom.v0   = v0;
1252:     geom.J    = J;
1253:     geom.detJ = detJ;
1254:     for (f = 0, v = 0; f < numFields; ++f) {
1255:       void * const ctx = ctxs ? ctxs[f] : NULL;

1257:       PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1258:       PetscFEGetNumComponents(fe, &numComp);
1259:       PetscFEGetDualSpace(fe, &sp);
1260:       PetscDualSpaceGetDimension(sp, &spDim);
1261:       for (d = 0; d < spDim; ++d) {
1262:         PetscDualSpaceApply(sp, d, geom, numComp, funcs[f], ctx, &values[v]);
1263:         v += numComp;
1264:       }
1265:     }
1266:     DMDAVecSetClosure(dm, section, localX, c, values, mode);
1267:   }
1268:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
1269:   PetscFree3(v0,J,detJ);
1270:   return(0);
1271: }

1275: /*@C
1276:   DMDAProjectFunction - This projects the given function into the function space provided.

1278:   Input Parameters:
1279: + dm      - The DM
1280: . funcs   - The coordinate functions to evaluate
1281: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
1282: - mode    - The insertion mode for values

1284:   Output Parameter:
1285: . X - vector

1287:   Level: developer

1289: .seealso: DMDAComputeL2Diff()
1290: @*/
1291: PetscErrorCode DMDAProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
1292: {
1293:   Vec            localX;

1298:   DMGetLocalVector(dm, &localX);
1299:   DMDAProjectFunctionLocal(dm, funcs, ctxs, mode, localX);
1300:   DMLocalToGlobalBegin(dm, localX, mode, X);
1301:   DMLocalToGlobalEnd(dm, localX, mode, X);
1302:   DMRestoreLocalVector(dm, &localX);
1303:   return(0);
1304: }

1308: /*@C
1309:   DMDAComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

1311:   Input Parameters:
1312: + dm    - The DM
1313: . funcs - The functions to evaluate for each field component
1314: . ctxs  - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
1315: - X     - The coefficient vector u_h

1317:   Output Parameter:
1318: . diff - The diff ||u - u_h||_2

1320:   Level: developer

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

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

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

1360:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1361:       void * const ctx = ctxs ? ctxs[field] : NULL;
1362:       const PetscReal *quadPoints, *quadWeights;
1363:       PetscReal       *basis;
1364:       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;

1366:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1367:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1368:       PetscFEGetDimension(fe, &numBasisFuncs);
1369:       PetscFEGetNumComponents(fe, &numBasisComps);
1370:       PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
1371:       if (debug) {
1372:         char title[1024];
1373:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1374:         DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
1375:       }
1376:       for (q = 0; q < numQuadPoints; ++q) {
1377:         for (d = 0; d < dim; d++) {
1378:           coords[d] = v0[d];
1379:           for (e = 0; e < dim; e++) {
1380:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1381:           }
1382:         }
1383:         (*funcs[field])(coords, funcVal, ctx);
1384:         for (fc = 0; fc < numBasisComps; ++fc) {
1385:           PetscScalar interpolant = 0.0;

1387:           for (f = 0; f < numBasisFuncs; ++f) {
1388:             const PetscInt fidx = f*numBasisComps+fc;
1389:             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
1390:           }
1391:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
1392:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
1393:         }
1394:       }
1395:       comp        += numBasisComps;
1396:       fieldOffset += numBasisFuncs*numBasisComps;
1397:     }
1398:     DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1399:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
1400:     localDiff += elemDiff;
1401:   }
1402:   PetscFree5(funcVal,coords,v0,J,invJ);
1403:   DMRestoreLocalVector(dm, &localX);
1404:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
1405:   *diff = PetscSqrtReal(*diff);
1406:   return(0);
1407: }

1411: /*@C
1412:   DMDAComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

1414:   Input Parameters:
1415: + dm    - The DM
1416: . funcs - The gradient functions to evaluate for each field component
1417: . ctxs  - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
1418: . X     - The coefficient vector u_h
1419: - n     - The vector to project along

1421:   Output Parameter:
1422: . diff - The diff ||(grad u - grad u_h) . n||_2

1424:   Level: developer

1426: .seealso: DMDAProjectFunction(), DMPlexComputeL2Diff()
1427: @*/
1428: PetscErrorCode DMDAComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
1429: {
1430:   const PetscInt  debug = 0;
1431:   PetscDS    prob;
1432:   PetscFE         fe;
1433:   PetscQuadrature quad;
1434:   PetscSection    section;
1435:   Vec             localX;
1436:   PetscScalar    *funcVal, *interpolantVec;
1437:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
1438:   PetscReal       localDiff = 0.0;
1439:   PetscInt        dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1440:   PetscErrorCode  ierr;

1443:   DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1444:   DMGetDS(dm, &prob);
1445:   PetscDSGetTotalComponents(prob, &Nc);
1446:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1447:   PetscFEGetQuadrature(fe, &quad);
1448:   DMGetDefaultSection(dm, &section);
1449:   PetscSectionGetNumFields(section, &numFields);
1450:   DMGetLocalVector(dm, &localX);
1451:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1452:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1453:   /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1454:   PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
1455:   DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1456:   for (c = cStart; c < cEnd; ++c) {
1457:     PetscScalar *x = NULL;
1458:     PetscReal    elemDiff = 0.0;

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

1464:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1465:       void * const ctx = ctxs ? ctxs[field] : NULL;
1466:       const PetscReal *quadPoints, *quadWeights;
1467:       PetscReal       *basisDer;
1468:       PetscInt         Nq, Nb, Nc, q, d, e, fc, f, g;

1470:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1471:       PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
1472:       PetscFEGetDimension(fe, &Nb);
1473:       PetscFEGetNumComponents(fe, &Nc);
1474:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1475:       if (debug) {
1476:         char title[1024];
1477:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1478:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1479:       }
1480:       for (q = 0; q < Nq; ++q) {
1481:         for (d = 0; d < dim; d++) {
1482:           coords[d] = v0[d];
1483:           for (e = 0; e < dim; e++) {
1484:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1485:           }
1486:         }
1487:         (*funcs[field])(coords, n, funcVal, ctx);
1488:         for (fc = 0; fc < Nc; ++fc) {
1489:           PetscScalar interpolant = 0.0;

1491:           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
1492:           for (f = 0; f < Nb; ++f) {
1493:             const PetscInt fidx = f*Nc+fc;

1495:             for (d = 0; d < dim; ++d) {
1496:               realSpaceDer[d] = 0.0;
1497:               for (g = 0; g < dim; ++g) {
1498:                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Nc+fidx)*dim+g];
1499:               }
1500:               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
1501:             }
1502:           }
1503:           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
1504:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
1505:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
1506:         }
1507:       }
1508:       comp        += Nc;
1509:       fieldOffset += Nb*Nc;
1510:     }
1511:     DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1512:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
1513:     localDiff += elemDiff;
1514:   }
1515:   PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
1516:   DMRestoreLocalVector(dm, &localX);
1517:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
1518:   *diff = PetscSqrtReal(*diff);
1519:   return(0);
1520: }

1522: /* ------------------------------------------------------------------- */

1526: /*@C
1527:      DMDAGetArray - Gets a work array for a DMDA

1529:     Input Parameter:
1530: +    da - information about my local patch
1531: -    ghosted - do you want arrays for the ghosted or nonghosted patch

1533:     Output Parameters:
1534: .    vptr - array data structured

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

1539:   Level: advanced

1541: .seealso: DMDARestoreArray()

1543: @*/
1544: PetscErrorCode  DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
1545: {
1547:   PetscInt       j,i,xs,ys,xm,ym,zs,zm;
1548:   char           *iarray_start;
1549:   void           **iptr = (void**)vptr;
1550:   DM_DA          *dd    = (DM_DA*)da->data;

1554:   if (ghosted) {
1555:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1556:       if (dd->arrayghostedin[i]) {
1557:         *iptr                 = dd->arrayghostedin[i];
1558:         iarray_start          = (char*)dd->startghostedin[i];
1559:         dd->arrayghostedin[i] = NULL;
1560:         dd->startghostedin[i] = NULL;

1562:         goto done;
1563:       }
1564:     }
1565:     xs = dd->Xs;
1566:     ys = dd->Ys;
1567:     zs = dd->Zs;
1568:     xm = dd->Xe-dd->Xs;
1569:     ym = dd->Ye-dd->Ys;
1570:     zm = dd->Ze-dd->Zs;
1571:   } else {
1572:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1573:       if (dd->arrayin[i]) {
1574:         *iptr          = dd->arrayin[i];
1575:         iarray_start   = (char*)dd->startin[i];
1576:         dd->arrayin[i] = NULL;
1577:         dd->startin[i] = NULL;

1579:         goto done;
1580:       }
1581:     }
1582:     xs = dd->xs;
1583:     ys = dd->ys;
1584:     zs = dd->zs;
1585:     xm = dd->xe-dd->xs;
1586:     ym = dd->ye-dd->ys;
1587:     zm = dd->ze-dd->zs;
1588:   }

1590:   switch (dd->dim) {
1591:   case 1: {
1592:     void *ptr;

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

1596:     ptr   = (void*)(iarray_start - xs*sizeof(PetscScalar));
1597:     *iptr = (void*)ptr;
1598:     break;
1599:   }
1600:   case 2: {
1601:     void **ptr;

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

1605:     ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
1606:     for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
1607:     *iptr = (void*)ptr;
1608:     break;
1609:   }
1610:   case 3: {
1611:     void ***ptr,**bptr;

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

1615:     ptr  = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1616:     bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1617:     for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1618:     for (i=zs; i<zs+zm; i++) {
1619:       for (j=ys; j<ys+ym; j++) {
1620:         ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1621:       }
1622:     }
1623:     *iptr = (void*)ptr;
1624:     break;
1625:   }
1626:   default:
1627:     SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
1628:   }

1630: done:
1631:   /* add arrays to the checked out list */
1632:   if (ghosted) {
1633:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1634:       if (!dd->arrayghostedout[i]) {
1635:         dd->arrayghostedout[i] = *iptr;
1636:         dd->startghostedout[i] = iarray_start;
1637:         break;
1638:       }
1639:     }
1640:   } else {
1641:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1642:       if (!dd->arrayout[i]) {
1643:         dd->arrayout[i] = *iptr;
1644:         dd->startout[i] = iarray_start;
1645:         break;
1646:       }
1647:     }
1648:   }
1649:   return(0);
1650: }

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

1657:     Input Parameter:
1658: +    da - information about my local patch
1659: .    ghosted - do you want arrays for the ghosted or nonghosted patch
1660: -    vptr - array data structured to be passed to ad_FormFunctionLocal()

1662:      Level: advanced

1664: .seealso: DMDAGetArray()

1666: @*/
1667: PetscErrorCode  DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1668: {
1669:   PetscInt i;
1670:   void     **iptr = (void**)vptr,*iarray_start = 0;
1671:   DM_DA    *dd    = (DM_DA*)da->data;

1675:   if (ghosted) {
1676:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1677:       if (dd->arrayghostedout[i] == *iptr) {
1678:         iarray_start           = dd->startghostedout[i];
1679:         dd->arrayghostedout[i] = NULL;
1680:         dd->startghostedout[i] = NULL;
1681:         break;
1682:       }
1683:     }
1684:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1685:       if (!dd->arrayghostedin[i]) {
1686:         dd->arrayghostedin[i] = *iptr;
1687:         dd->startghostedin[i] = iarray_start;
1688:         break;
1689:       }
1690:     }
1691:   } else {
1692:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1693:       if (dd->arrayout[i] == *iptr) {
1694:         iarray_start    = dd->startout[i];
1695:         dd->arrayout[i] = NULL;
1696:         dd->startout[i] = NULL;
1697:         break;
1698:       }
1699:     }
1700:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1701:       if (!dd->arrayin[i]) {
1702:         dd->arrayin[i] = *iptr;
1703:         dd->startin[i] = iarray_start;
1704:         break;
1705:       }
1706:     }
1707:   }
1708:   return(0);
1709: }