Actual source code: dalocal.c

petsc-master 2014-10-19
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  && da->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 = dm->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:   const PetscInt dim = dm->dim;
142:   DMDALocalInfo  info;

148:   DMDAGetLocalInfo(dm, &info);
149:   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);}
150:   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);}
151:   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);}
152:   *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0);
153:   return(0);
154: }

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

169:   if (numVerticesX) {
171:     *numVerticesX = nVx;
172:   }
173:   if (numVerticesY) {
175:     *numVerticesY = nVy;
176:   }
177:   if (numVerticesZ) {
179:     *numVerticesZ = nVz;
180:   }
181:   if (numVertices) {
183:     *numVertices = nV;
184:   }
185:   return(0);
186: }

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

203:   if (numXFacesX) {
205:     *numXFacesX = nxF;
206:   }
207:   if (numXFaces) {
209:     *numXFaces = nXF;
210:   }
211:   if (numYFacesY) {
213:     *numYFacesY = nyF;
214:   }
215:   if (numYFaces) {
217:     *numYFaces = nYF;
218:   }
219:   if (numZFacesZ) {
221:     *numZFacesZ = nzF;
222:   }
223:   if (numZFaces) {
225:     *numZFaces = nZF;
226:   }
227:   return(0);
228: }

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

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

266: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
267: {
268:   const PetscInt dim = dm->dim;
269:   PetscInt       nC, nV, nXF, nYF, nZF;

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

300: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
301: {
302:   const PetscInt dim = dm->dim;
303:   PetscInt       nC, nV, nXF, nYF, nZF;

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

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

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

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

360:         (*cone)[0] = fy*nVx + fx + nC;
361:         (*cone)[1] = fy*nVx + fx + 1 + nC;
362:       } else if (p < nC+nV+nXF+nYF) {
363:         const PetscInt fx = (p - nC+nV+nXF) / nyF;
364:         const PetscInt fy = (p - nC+nV+nXF) % nyF;

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

380: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
381: {

385:   DMGetWorkArray(dm, 6, PETSC_INT, cone);
386:   return(0);
387: }

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

395:   Input Parameters:
396: + dm- The DMDA
397: . numFields - The number of fields
398: . numComp - The number of components in each field
399: . numDof - The number of dofs per dimension for each field
400: . numFaceDof - The number of dofs per face for each field and direction, or NULL

402:   Level: developer

404:   Note:
405:   The default DMDA numbering is as follows:

407:     - Cells:    [0,             nC)
408:     - Vertices: [nC,            nC+nV)
409:     - X-Faces:  [nC+nV,         nC+nV+nXF)         normal is +- x-dir
410:     - Y-Faces:  [nC+nV+nXF,     nC+nV+nXF+nYF)     normal is +- y-dir
411:     - Z-Faces:  [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir

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

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

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

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

720:         if (neighbor >= 0 && neighbor < rank) {
721:           if (xp < 0) { /* left */
722:             if (yp < 0) { /* bottom */
723:               if (zp < 0) { /* back */
724:                 const PetscInt localVertex  = (      0*nVy +     0)*nVx +     0 + nC; /* left bottom back vertex */
725:                 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

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

737:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
738:                   localPoints[nL]        = localVertex;
739:                   remotePoints[nL].rank  = neighbor;
740:                   remotePoints[nL].index = remoteVertex;
741:                   ++nL;
742:                 }
743:               } else {
744:                 for (zv = 0; zv < nVz; ++zv) {
745:                   const PetscInt localVertex  = (zv*nVy +     0)*nVx +     0 + nC; /* left bottom vertices */
746:                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*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:                 }
755:               }
756:             } else if (yp > 0) { /* top */
757:               if (zp < 0) { /* back */
758:                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx +     0 + nC; /* left top back vertex */
759:                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

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

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

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

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

808:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
809:                     localPoints[nL]        = localVertex;
810:                     remotePoints[nL].rank  = neighbor;
811:                     remotePoints[nL].index = remoteVertex;
812:                     ++nL;
813:                   }
814:                 }
815:               } else {
816:                 for (zv = 0; zv < nVz; ++zv) {
817:                   for (yv = 0; yv < nVy; ++yv) {
818:                     const PetscInt localVertex  = (zv*nVy + yv)*nVx +     0 + nC; /* left vertices */
819:                     const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */

821:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
822:                       localPoints[nL]        = localVertex;
823:                       remotePoints[nL].rank  = neighbor;
824:                       remotePoints[nL].index = remoteVertex;
825:                       ++nL;
826:                     }
827:                   }
828:                 }
829: #if 0
830:                 for (xf = 0; xf < nxF; ++xf) {
831:                   /* THIS IS WRONG */
832:                   const PetscInt localFace  = 0 + nC+nV; /* left faces */
833:                   const PetscInt remoteFace = 0 + nC+nV;

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

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

860:                 if (!PetscBTLookupSet(isLeaf, localVertex)) {
861:                   localPoints[nL]        = localVertex;
862:                   remotePoints[nL].rank  = neighbor;
863:                   remotePoints[nL].index = remoteVertex;
864:                   ++nL;
865:                 }
866:               } else {
867:                 nleavesCheck += nVz;
868:                 for (zv = 0; zv < nVz; ++zv) {
869:                   const PetscInt localVertex  = (zv*nVy +     0)*nVx + nVx-1 + nC; /* right bottom vertices */
870:                   const PetscInt remoteVertex = (zv*nVy + nVy-1)*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:                 }
879:               }
880:             } else if (yp > 0) { /* top */
881:               if (zp < 0) { /* back */
882:                 const PetscInt localVertex  = (      0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
883:                 const PetscInt remoteVertex = ((nVz-1)*nVy +     0)*nVx +     0 + nC; /* TODO: Correct this for neighbor sizes */

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

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

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

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

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

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

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

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

988:                   if (!PetscBTLookupSet(isLeaf, localVertex)) {
989:                     localPoints[nL]        = localVertex;
990:                     remotePoints[nL].rank  = neighbor;
991:                     remotePoints[nL].index = remoteVertex;
992:                     ++nL;
993:                   }
994:                 }
995:               } else {
996:                 for (zv = 0; zv < nVz; ++zv) {
997:                   for (xv = 0; xv < nVx; ++xv) {
998:                     const PetscInt localVertex  = (zv*nVy +     0)*nVx + xv + nC; /* bottom vertices */
999:                     const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */

1001:                     if (!PetscBTLookupSet(isLeaf, localVertex)) {
1002:                       localPoints[nL]        = localVertex;
1003:                       remotePoints[nL].rank  = neighbor;
1004:                       remotePoints[nL].index = remoteVertex;
1005:                       ++nL;
1006:                     }
1007:                   }
1008:                 }
1009: #if 0
1010:                 for (yf = 0; yf < nyF; ++yf) {
1011:                   /* THIS IS WRONG */
1012:                   const PetscInt localFace  = 0 + nC+nV; /* bottom faces */
1013:                   const PetscInt remoteFace = 0 + nC+nV;

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

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

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

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

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

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

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

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

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

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

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

1189:     ind[0] = 0;
1190:     ind[1] = 0;
1191:     ind[2] = k + da->zs;
1192:     for (j = 0; j < nVy; ++j) {
1193:       ind[1] = j + da->ys;
1194:       for (i = 0; i < nVx; ++i) {
1195:         const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;

1197:         PetscSectionGetOffset(section, vertex, &off);
1198:         ind[0] = i + da->xs;
1199:         for (d = 0; d < dim; ++d) {
1200:           coords[off+d] = h[d]*ind[d];
1201:         }
1202:       }
1203:     }
1204:   }
1205:   VecRestoreArray(coordinates, &coords);
1206:   DMSetCoordinateSection(dm, PETSC_DETERMINE, section);
1207:   DMSetCoordinatesLocal(dm, coordinates);
1208:   PetscSectionDestroy(&section);
1209:   VecDestroy(&coordinates);
1210:   return(0);
1211: }

1215: PetscErrorCode DMDAProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1216: {
1217:   PetscDS    prob;
1218:   PetscFE         fe;
1219:   PetscDualSpace  sp;
1220:   PetscQuadrature q;
1221:   PetscSection    section;
1222:   PetscScalar    *values;
1223:   PetscInt        numFields, numComp, numPoints, dim, spDim, totDim, numValues, cStart, cEnd, f, c, v, d;
1224:   PetscErrorCode  ierr;

1227:   DMGetDS(dm, &prob);
1228:   PetscDSGetTotalDimension(prob, &totDim);
1229:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1230:   PetscFEGetQuadrature(fe, &q);
1231:   DMGetDefaultSection(dm, &section);
1232:   PetscSectionGetNumFields(section, &numFields);
1233:   DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1234:   DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1235:   DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
1236:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
1237:   DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
1238:   PetscQuadratureGetData(q, NULL, &numPoints, NULL, NULL);
1239:   for (c = cStart; c < cEnd; ++c) {
1240:     PetscFECellGeom geom;

1242:     DMDAComputeCellGeometryFEM(dm, c, q, geom.v0, geom.J, NULL, &geom.detJ);
1243:     for (f = 0, v = 0; f < numFields; ++f) {
1244:       void * const ctx = ctxs ? ctxs[f] : NULL;

1246:       PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1247:       PetscFEGetNumComponents(fe, &numComp);
1248:       PetscFEGetDualSpace(fe, &sp);
1249:       PetscDualSpaceGetDimension(sp, &spDim);
1250:       for (d = 0; d < spDim; ++d) {
1251:         PetscDualSpaceApply(sp, d, &geom, numComp, funcs[f], ctx, &values[v]);
1252:         v += numComp;
1253:       }
1254:     }
1255:     DMDAVecSetClosure(dm, section, localX, c, values, mode);
1256:   }
1257:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
1258:   return(0);
1259: }

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

1266:   Input Parameters:
1267: + dm      - The DM
1268: . funcs   - The coordinate functions to evaluate
1269: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
1270: - mode    - The insertion mode for values

1272:   Output Parameter:
1273: . X - vector

1275:   Level: developer

1277: .seealso: DMDAComputeL2Diff()
1278: @*/
1279: PetscErrorCode DMDAProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
1280: {
1281:   Vec            localX;

1286:   DMGetLocalVector(dm, &localX);
1287:   DMDAProjectFunctionLocal(dm, funcs, ctxs, mode, localX);
1288:   DMLocalToGlobalBegin(dm, localX, mode, X);
1289:   DMLocalToGlobalEnd(dm, localX, mode, X);
1290:   DMRestoreLocalVector(dm, &localX);
1291:   return(0);
1292: }

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

1299:   Input Parameters:
1300: + dm    - The DM
1301: . funcs - The functions to evaluate for each field component
1302: . ctxs  - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
1303: - X     - The coefficient vector u_h

1305:   Output Parameter:
1306: . diff - The diff ||u - u_h||_2

1308:   Level: developer

1310: .seealso: DMDAProjectFunction(), DMDAComputeL2GradientDiff()
1311: @*/
1312: PetscErrorCode DMDAComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1313: {
1314:   const PetscInt  debug = 0;
1315:   PetscDS    prob;
1316:   PetscFE         fe;
1317:   PetscQuadrature quad;
1318:   PetscSection    section;
1319:   Vec             localX;
1320:   PetscScalar    *funcVal;
1321:   PetscReal      *coords, *v0, *J, *invJ, detJ;
1322:   PetscReal       localDiff = 0.0;
1323:   PetscInt        dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1324:   PetscErrorCode  ierr;

1327:   DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1328:   DMGetDS(dm, &prob);
1329:   PetscDSGetTotalComponents(prob, &Nc);
1330:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1331:   PetscFEGetQuadrature(fe, &quad);
1332:   DMGetDefaultSection(dm, &section);
1333:   PetscSectionGetNumFields(section, &numFields);
1334:   DMGetLocalVector(dm, &localX);
1335:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1336:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1337:   /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1338:   PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1339:   DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1340:   for (c = cStart; c < cEnd; ++c) {
1341:     PetscScalar *x = NULL;
1342:     PetscReal    elemDiff = 0.0;

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

1348:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1349:       void * const ctx = ctxs ? ctxs[field] : NULL;
1350:       const PetscReal *quadPoints, *quadWeights;
1351:       PetscReal       *basis;
1352:       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;

1354:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1355:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1356:       PetscFEGetDimension(fe, &numBasisFuncs);
1357:       PetscFEGetNumComponents(fe, &numBasisComps);
1358:       PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
1359:       if (debug) {
1360:         char title[1024];
1361:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1362:         DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
1363:       }
1364:       for (q = 0; q < numQuadPoints; ++q) {
1365:         for (d = 0; d < dim; d++) {
1366:           coords[d] = v0[d];
1367:           for (e = 0; e < dim; e++) {
1368:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1369:           }
1370:         }
1371:         (*funcs[field])(coords, funcVal, ctx);
1372:         for (fc = 0; fc < numBasisComps; ++fc) {
1373:           PetscScalar interpolant = 0.0;

1375:           for (f = 0; f < numBasisFuncs; ++f) {
1376:             const PetscInt fidx = f*numBasisComps+fc;
1377:             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
1378:           }
1379:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
1380:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
1381:         }
1382:       }
1383:       comp        += numBasisComps;
1384:       fieldOffset += numBasisFuncs*numBasisComps;
1385:     }
1386:     DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1387:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
1388:     localDiff += elemDiff;
1389:   }
1390:   PetscFree5(funcVal,coords,v0,J,invJ);
1391:   DMRestoreLocalVector(dm, &localX);
1392:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
1393:   *diff = PetscSqrtReal(*diff);
1394:   return(0);
1395: }

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

1402:   Input Parameters:
1403: + dm    - The DM
1404: . funcs - The gradient functions to evaluate for each field component
1405: . ctxs  - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
1406: . X     - The coefficient vector u_h
1407: - n     - The vector to project along

1409:   Output Parameter:
1410: . diff - The diff ||(grad u - grad u_h) . n||_2

1412:   Level: developer

1414: .seealso: DMDAProjectFunction(), DMPlexComputeL2Diff()
1415: @*/
1416: PetscErrorCode DMDAComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
1417: {
1418:   const PetscInt  debug = 0;
1419:   PetscDS    prob;
1420:   PetscFE         fe;
1421:   PetscQuadrature quad;
1422:   PetscSection    section;
1423:   Vec             localX;
1424:   PetscScalar    *funcVal, *interpolantVec;
1425:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
1426:   PetscReal       localDiff = 0.0;
1427:   PetscInt        dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1428:   PetscErrorCode  ierr;

1431:   DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1432:   DMGetDS(dm, &prob);
1433:   PetscDSGetTotalComponents(prob, &Nc);
1434:   PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1435:   PetscFEGetQuadrature(fe, &quad);
1436:   DMGetDefaultSection(dm, &section);
1437:   PetscSectionGetNumFields(section, &numFields);
1438:   DMGetLocalVector(dm, &localX);
1439:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1440:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1441:   /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1442:   PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
1443:   DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1444:   for (c = cStart; c < cEnd; ++c) {
1445:     PetscScalar *x = NULL;
1446:     PetscReal    elemDiff = 0.0;

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

1452:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1453:       void * const ctx = ctxs ? ctxs[field] : NULL;
1454:       const PetscReal *quadPoints, *quadWeights;
1455:       PetscReal       *basisDer;
1456:       PetscInt         Nq, Nb, Nc, q, d, e, fc, f, g;

1458:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1459:       PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
1460:       PetscFEGetDimension(fe, &Nb);
1461:       PetscFEGetNumComponents(fe, &Nc);
1462:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1463:       if (debug) {
1464:         char title[1024];
1465:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1466:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1467:       }
1468:       for (q = 0; q < Nq; ++q) {
1469:         for (d = 0; d < dim; d++) {
1470:           coords[d] = v0[d];
1471:           for (e = 0; e < dim; e++) {
1472:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1473:           }
1474:         }
1475:         (*funcs[field])(coords, n, funcVal, ctx);
1476:         for (fc = 0; fc < Nc; ++fc) {
1477:           PetscScalar interpolant = 0.0;

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

1483:             for (d = 0; d < dim; ++d) {
1484:               realSpaceDer[d] = 0.0;
1485:               for (g = 0; g < dim; ++g) {
1486:                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Nc+fidx)*dim+g];
1487:               }
1488:               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
1489:             }
1490:           }
1491:           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
1492:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
1493:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
1494:         }
1495:       }
1496:       comp        += Nc;
1497:       fieldOffset += Nb*Nc;
1498:     }
1499:     DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1500:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
1501:     localDiff += elemDiff;
1502:   }
1503:   PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
1504:   DMRestoreLocalVector(dm, &localX);
1505:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
1506:   *diff = PetscSqrtReal(*diff);
1507:   return(0);
1508: }

1510: /* ------------------------------------------------------------------- */

1514: /*@C
1515:      DMDAGetArray - Gets a work array for a DMDA

1517:     Input Parameter:
1518: +    da - information about my local patch
1519: -    ghosted - do you want arrays for the ghosted or nonghosted patch

1521:     Output Parameters:
1522: .    vptr - array data structured

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

1527:   Level: advanced

1529: .seealso: DMDARestoreArray()

1531: @*/
1532: PetscErrorCode  DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
1533: {
1535:   PetscInt       j,i,xs,ys,xm,ym,zs,zm;
1536:   char           *iarray_start;
1537:   void           **iptr = (void**)vptr;
1538:   DM_DA          *dd    = (DM_DA*)da->data;

1542:   if (ghosted) {
1543:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1544:       if (dd->arrayghostedin[i]) {
1545:         *iptr                 = dd->arrayghostedin[i];
1546:         iarray_start          = (char*)dd->startghostedin[i];
1547:         dd->arrayghostedin[i] = NULL;
1548:         dd->startghostedin[i] = NULL;

1550:         goto done;
1551:       }
1552:     }
1553:     xs = dd->Xs;
1554:     ys = dd->Ys;
1555:     zs = dd->Zs;
1556:     xm = dd->Xe-dd->Xs;
1557:     ym = dd->Ye-dd->Ys;
1558:     zm = dd->Ze-dd->Zs;
1559:   } else {
1560:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1561:       if (dd->arrayin[i]) {
1562:         *iptr          = dd->arrayin[i];
1563:         iarray_start   = (char*)dd->startin[i];
1564:         dd->arrayin[i] = NULL;
1565:         dd->startin[i] = NULL;

1567:         goto done;
1568:       }
1569:     }
1570:     xs = dd->xs;
1571:     ys = dd->ys;
1572:     zs = dd->zs;
1573:     xm = dd->xe-dd->xs;
1574:     ym = dd->ye-dd->ys;
1575:     zm = dd->ze-dd->zs;
1576:   }

1578:   switch (da->dim) {
1579:   case 1: {
1580:     void *ptr;

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

1584:     ptr   = (void*)(iarray_start - xs*sizeof(PetscScalar));
1585:     *iptr = (void*)ptr;
1586:     break;
1587:   }
1588:   case 2: {
1589:     void **ptr;

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

1593:     ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
1594:     for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
1595:     *iptr = (void*)ptr;
1596:     break;
1597:   }
1598:   case 3: {
1599:     void ***ptr,**bptr;

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

1603:     ptr  = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1604:     bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1605:     for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1606:     for (i=zs; i<zs+zm; i++) {
1607:       for (j=ys; j<ys+ym; j++) {
1608:         ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1609:       }
1610:     }
1611:     *iptr = (void*)ptr;
1612:     break;
1613:   }
1614:   default:
1615:     SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
1616:   }

1618: done:
1619:   /* add arrays to the checked out list */
1620:   if (ghosted) {
1621:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1622:       if (!dd->arrayghostedout[i]) {
1623:         dd->arrayghostedout[i] = *iptr;
1624:         dd->startghostedout[i] = iarray_start;
1625:         break;
1626:       }
1627:     }
1628:   } else {
1629:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1630:       if (!dd->arrayout[i]) {
1631:         dd->arrayout[i] = *iptr;
1632:         dd->startout[i] = iarray_start;
1633:         break;
1634:       }
1635:     }
1636:   }
1637:   return(0);
1638: }

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

1645:     Input Parameter:
1646: +    da - information about my local patch
1647: .    ghosted - do you want arrays for the ghosted or nonghosted patch
1648: -    vptr - array data structured to be passed to ad_FormFunctionLocal()

1650:      Level: advanced

1652: .seealso: DMDAGetArray()

1654: @*/
1655: PetscErrorCode  DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1656: {
1657:   PetscInt i;
1658:   void     **iptr = (void**)vptr,*iarray_start = 0;
1659:   DM_DA    *dd    = (DM_DA*)da->data;

1663:   if (ghosted) {
1664:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1665:       if (dd->arrayghostedout[i] == *iptr) {
1666:         iarray_start           = dd->startghostedout[i];
1667:         dd->arrayghostedout[i] = NULL;
1668:         dd->startghostedout[i] = NULL;
1669:         break;
1670:       }
1671:     }
1672:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1673:       if (!dd->arrayghostedin[i]) {
1674:         dd->arrayghostedin[i] = *iptr;
1675:         dd->startghostedin[i] = iarray_start;
1676:         break;
1677:       }
1678:     }
1679:   } else {
1680:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1681:       if (dd->arrayout[i] == *iptr) {
1682:         iarray_start    = dd->startout[i];
1683:         dd->arrayout[i] = NULL;
1684:         dd->startout[i] = NULL;
1685:         break;
1686:       }
1687:     }
1688:     for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1689:       if (!dd->arrayin[i]) {
1690:         dd->arrayin[i] = *iptr;
1691:         dd->startin[i] = iarray_start;
1692:         break;
1693:       }
1694:     }
1695:   }
1696:   return(0);
1697: }