Actual source code: dalocal.c
petsc-3.7.5 2017-01-01
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), §ion);
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: PetscMalloc1(nleaves,&localPoints);
713: PetscMalloc1(nleaves,&remotePoints);
714: for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
715: for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
716: for (xn = 0; xn < 3; ++xn) {
717: const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
718: const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
719: PetscInt xv, yv, zv;
721: if (neighbor >= 0 && neighbor < rank) {
722: if (xp < 0) { /* left */
723: if (yp < 0) { /* bottom */
724: if (zp < 0) { /* back */
725: const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* left bottom back vertex */
726: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
728: if (!PetscBTLookupSet(isLeaf, localVertex)) {
729: localPoints[nL] = localVertex;
730: remotePoints[nL].rank = neighbor;
731: remotePoints[nL].index = remoteVertex;
732: ++nL;
733: }
734: } else if (zp > 0) { /* front */
735: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* left bottom front vertex */
736: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
738: if (!PetscBTLookupSet(isLeaf, localVertex)) {
739: localPoints[nL] = localVertex;
740: remotePoints[nL].rank = neighbor;
741: remotePoints[nL].index = remoteVertex;
742: ++nL;
743: }
744: } else {
745: for (zv = 0; zv < nVz; ++zv) {
746: const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; /* left bottom vertices */
747: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
749: if (!PetscBTLookupSet(isLeaf, localVertex)) {
750: localPoints[nL] = localVertex;
751: remotePoints[nL].rank = neighbor;
752: remotePoints[nL].index = remoteVertex;
753: ++nL;
754: }
755: }
756: }
757: } else if (yp > 0) { /* top */
758: if (zp < 0) { /* back */
759: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* left top back vertex */
760: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
762: if (!PetscBTLookupSet(isLeaf, localVertex)) {
763: localPoints[nL] = localVertex;
764: remotePoints[nL].rank = neighbor;
765: remotePoints[nL].index = remoteVertex;
766: ++nL;
767: }
768: } else if (zp > 0) { /* front */
769: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* left top front vertex */
770: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
772: if (!PetscBTLookupSet(isLeaf, localVertex)) {
773: localPoints[nL] = localVertex;
774: remotePoints[nL].rank = neighbor;
775: remotePoints[nL].index = remoteVertex;
776: ++nL;
777: }
778: } else {
779: for (zv = 0; zv < nVz; ++zv) {
780: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* left top vertices */
781: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
783: if (!PetscBTLookupSet(isLeaf, localVertex)) {
784: localPoints[nL] = localVertex;
785: remotePoints[nL].rank = neighbor;
786: remotePoints[nL].index = remoteVertex;
787: ++nL;
788: }
789: }
790: }
791: } else {
792: if (zp < 0) { /* back */
793: for (yv = 0; yv < nVy; ++yv) {
794: const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* left back vertices */
795: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
797: if (!PetscBTLookupSet(isLeaf, localVertex)) {
798: localPoints[nL] = localVertex;
799: remotePoints[nL].rank = neighbor;
800: remotePoints[nL].index = remoteVertex;
801: ++nL;
802: }
803: }
804: } else if (zp > 0) { /* front */
805: for (yv = 0; yv < nVy; ++yv) {
806: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* left front vertices */
807: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
809: if (!PetscBTLookupSet(isLeaf, localVertex)) {
810: localPoints[nL] = localVertex;
811: remotePoints[nL].rank = neighbor;
812: remotePoints[nL].index = remoteVertex;
813: ++nL;
814: }
815: }
816: } else {
817: for (zv = 0; zv < nVz; ++zv) {
818: for (yv = 0; yv < nVy; ++yv) {
819: const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; /* left vertices */
820: const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
822: if (!PetscBTLookupSet(isLeaf, localVertex)) {
823: localPoints[nL] = localVertex;
824: remotePoints[nL].rank = neighbor;
825: remotePoints[nL].index = remoteVertex;
826: ++nL;
827: }
828: }
829: }
830: #if 0
831: for (xf = 0; xf < nxF; ++xf) {
832: /* THIS IS WRONG */
833: const PetscInt localFace = 0 + nC+nV; /* left faces */
834: const PetscInt remoteFace = 0 + nC+nV;
836: if (!PetscBTLookupSet(isLeaf, localFace)) {
837: localPoints[nL] = localFace;
838: remotePoints[nL].rank = neighbor;
839: remotePoints[nL].index = remoteFace;
840: }
841: }
842: #endif
843: }
844: }
845: } else if (xp > 0) { /* right */
846: if (yp < 0) { /* bottom */
847: if (zp < 0) { /* back */
848: const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* right bottom back vertex */
849: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
851: if (!PetscBTLookupSet(isLeaf, localVertex)) {
852: localPoints[nL] = localVertex;
853: remotePoints[nL].rank = neighbor;
854: remotePoints[nL].index = remoteVertex;
855: ++nL;
856: }
857: } else if (zp > 0) { /* front */
858: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* right bottom front vertex */
859: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
861: if (!PetscBTLookupSet(isLeaf, localVertex)) {
862: localPoints[nL] = localVertex;
863: remotePoints[nL].rank = neighbor;
864: remotePoints[nL].index = remoteVertex;
865: ++nL;
866: }
867: } else {
868: nleavesCheck += nVz;
869: for (zv = 0; zv < nVz; ++zv) {
870: const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* right bottom vertices */
871: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
873: if (!PetscBTLookupSet(isLeaf, localVertex)) {
874: localPoints[nL] = localVertex;
875: remotePoints[nL].rank = neighbor;
876: remotePoints[nL].index = remoteVertex;
877: ++nL;
878: }
879: }
880: }
881: } else if (yp > 0) { /* top */
882: if (zp < 0) { /* back */
883: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
884: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
886: if (!PetscBTLookupSet(isLeaf, localVertex)) {
887: localPoints[nL] = localVertex;
888: remotePoints[nL].rank = neighbor;
889: remotePoints[nL].index = remoteVertex;
890: ++nL;
891: }
892: } else if (zp > 0) { /* front */
893: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
894: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
896: if (!PetscBTLookupSet(isLeaf, localVertex)) {
897: localPoints[nL] = localVertex;
898: remotePoints[nL].rank = neighbor;
899: remotePoints[nL].index = remoteVertex;
900: ++nL;
901: }
902: } else {
903: for (zv = 0; zv < nVz; ++zv) {
904: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
905: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
907: if (!PetscBTLookupSet(isLeaf, localVertex)) {
908: localPoints[nL] = localVertex;
909: remotePoints[nL].rank = neighbor;
910: remotePoints[nL].index = remoteVertex;
911: ++nL;
912: }
913: }
914: }
915: } else {
916: if (zp < 0) { /* back */
917: for (yv = 0; yv < nVy; ++yv) {
918: const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
919: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
921: if (!PetscBTLookupSet(isLeaf, localVertex)) {
922: localPoints[nL] = localVertex;
923: remotePoints[nL].rank = neighbor;
924: remotePoints[nL].index = remoteVertex;
925: ++nL;
926: }
927: }
928: } else if (zp > 0) { /* front */
929: for (yv = 0; yv < nVy; ++yv) {
930: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
931: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
933: if (!PetscBTLookupSet(isLeaf, localVertex)) {
934: localPoints[nL] = localVertex;
935: remotePoints[nL].rank = neighbor;
936: remotePoints[nL].index = remoteVertex;
937: ++nL;
938: }
939: }
940: } else {
941: for (zv = 0; zv < nVz; ++zv) {
942: for (yv = 0; yv < nVy; ++yv) {
943: const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
944: const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
946: if (!PetscBTLookupSet(isLeaf, localVertex)) {
947: localPoints[nL] = localVertex;
948: remotePoints[nL].rank = neighbor;
949: remotePoints[nL].index = remoteVertex;
950: ++nL;
951: }
952: }
953: }
954: #if 0
955: for (xf = 0; xf < nxF; ++xf) {
956: /* THIS IS WRONG */
957: const PetscInt localFace = 0 + nC+nV; /* right faces */
958: const PetscInt remoteFace = 0 + nC+nV;
960: if (!PetscBTLookupSet(isLeaf, localFace)) {
961: localPoints[nL] = localFace;
962: remotePoints[nL].rank = neighbor;
963: remotePoints[nL].index = remoteFace;
964: ++nL;
965: }
966: }
967: #endif
968: }
969: }
970: } else {
971: if (yp < 0) { /* bottom */
972: if (zp < 0) { /* back */
973: for (xv = 0; xv < nVx; ++xv) {
974: const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; /* bottom back vertices */
975: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
977: if (!PetscBTLookupSet(isLeaf, localVertex)) {
978: localPoints[nL] = localVertex;
979: remotePoints[nL].rank = neighbor;
980: remotePoints[nL].index = remoteVertex;
981: ++nL;
982: }
983: }
984: } else if (zp > 0) { /* front */
985: for (xv = 0; xv < nVx; ++xv) {
986: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* bottom front vertices */
987: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
989: if (!PetscBTLookupSet(isLeaf, localVertex)) {
990: localPoints[nL] = localVertex;
991: remotePoints[nL].rank = neighbor;
992: remotePoints[nL].index = remoteVertex;
993: ++nL;
994: }
995: }
996: } else {
997: for (zv = 0; zv < nVz; ++zv) {
998: for (xv = 0; xv < nVx; ++xv) {
999: const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; /* bottom vertices */
1000: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1002: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1003: localPoints[nL] = localVertex;
1004: remotePoints[nL].rank = neighbor;
1005: remotePoints[nL].index = remoteVertex;
1006: ++nL;
1007: }
1008: }
1009: }
1010: #if 0
1011: for (yf = 0; yf < nyF; ++yf) {
1012: /* THIS IS WRONG */
1013: const PetscInt localFace = 0 + nC+nV; /* bottom faces */
1014: const PetscInt remoteFace = 0 + nC+nV;
1016: if (!PetscBTLookupSet(isLeaf, localFace)) {
1017: localPoints[nL] = localFace;
1018: remotePoints[nL].rank = neighbor;
1019: remotePoints[nL].index = remoteFace;
1020: ++nL;
1021: }
1022: }
1023: #endif
1024: }
1025: } else if (yp > 0) { /* top */
1026: if (zp < 0) { /* back */
1027: for (xv = 0; xv < nVx; ++xv) {
1028: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
1029: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1031: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1032: localPoints[nL] = localVertex;
1033: remotePoints[nL].rank = neighbor;
1034: remotePoints[nL].index = remoteVertex;
1035: ++nL;
1036: }
1037: }
1038: } else if (zp > 0) { /* front */
1039: for (xv = 0; xv < nVx; ++xv) {
1040: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
1041: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1043: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1044: localPoints[nL] = localVertex;
1045: remotePoints[nL].rank = neighbor;
1046: remotePoints[nL].index = remoteVertex;
1047: ++nL;
1048: }
1049: }
1050: } else {
1051: for (zv = 0; zv < nVz; ++zv) {
1052: for (xv = 0; xv < nVx; ++xv) {
1053: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
1054: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1056: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1057: localPoints[nL] = localVertex;
1058: remotePoints[nL].rank = neighbor;
1059: remotePoints[nL].index = remoteVertex;
1060: ++nL;
1061: }
1062: }
1063: }
1064: #if 0
1065: for (yf = 0; yf < nyF; ++yf) {
1066: /* THIS IS WRONG */
1067: const PetscInt localFace = 0 + nC+nV; /* top faces */
1068: const PetscInt remoteFace = 0 + nC+nV;
1070: if (!PetscBTLookupSet(isLeaf, localFace)) {
1071: localPoints[nL] = localFace;
1072: remotePoints[nL].rank = neighbor;
1073: remotePoints[nL].index = remoteFace;
1074: ++nL;
1075: }
1076: }
1077: #endif
1078: }
1079: } else {
1080: if (zp < 0) { /* back */
1081: for (yv = 0; yv < nVy; ++yv) {
1082: for (xv = 0; xv < nVx; ++xv) {
1083: const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; /* back vertices */
1084: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1086: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1087: localPoints[nL] = localVertex;
1088: remotePoints[nL].rank = neighbor;
1089: remotePoints[nL].index = remoteVertex;
1090: ++nL;
1091: }
1092: }
1093: }
1094: #if 0
1095: for (zf = 0; zf < nzF; ++zf) {
1096: /* THIS IS WRONG */
1097: const PetscInt localFace = 0 + nC+nV; /* back faces */
1098: const PetscInt remoteFace = 0 + nC+nV;
1100: if (!PetscBTLookupSet(isLeaf, localFace)) {
1101: localPoints[nL] = localFace;
1102: remotePoints[nL].rank = neighbor;
1103: remotePoints[nL].index = remoteFace;
1104: ++nL;
1105: }
1106: }
1107: #endif
1108: } else if (zp > 0) { /* front */
1109: for (yv = 0; yv < nVy; ++yv) {
1110: for (xv = 0; xv < nVx; ++xv) {
1111: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
1112: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1114: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1115: localPoints[nL] = localVertex;
1116: remotePoints[nL].rank = neighbor;
1117: remotePoints[nL].index = remoteVertex;
1118: ++nL;
1119: }
1120: }
1121: }
1122: #if 0
1123: for (zf = 0; zf < nzF; ++zf) {
1124: /* THIS IS WRONG */
1125: const PetscInt localFace = 0 + nC+nV; /* front faces */
1126: const PetscInt remoteFace = 0 + nC+nV;
1128: if (!PetscBTLookupSet(isLeaf, localFace)) {
1129: localPoints[nL] = localFace;
1130: remotePoints[nL].rank = neighbor;
1131: remotePoints[nL].index = remoteFace;
1132: ++nL;
1133: }
1134: }
1135: #endif
1136: } else {
1137: /* Nothing is shared from the interior */
1138: }
1139: }
1140: }
1141: }
1142: }
1143: }
1144: }
1145: PetscBTDestroy(&isLeaf);
1146: /* Remove duplication in leaf determination */
1147: 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);
1148: PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);
1149: PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1150: DMSetPointSF(dm, sf);
1151: PetscSFDestroy(&sf);
1152: *s = section;
1153: return(0);
1154: }
1158: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
1159: {
1160: DM_DA *da = (DM_DA *) dm->data;
1161: Vec coordinates;
1162: PetscSection section;
1163: PetscScalar *coords;
1164: PetscReal h[3];
1165: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
1170: DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);
1171: h[0] = (xu - xl)/M;
1172: h[1] = (yu - yl)/N;
1173: h[2] = (zu - zl)/P;
1174: DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
1175: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
1176: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);
1177: PetscSectionSetNumFields(section, 1);
1178: PetscSectionSetFieldComponents(section, 0, dim);
1179: PetscSectionSetChart(section, vStart, vEnd);
1180: for (v = vStart; v < vEnd; ++v) {
1181: PetscSectionSetDof(section, v, dim);
1182: }
1183: PetscSectionSetUp(section);
1184: PetscSectionGetStorageSize(section, &size);
1185: VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
1186: PetscObjectSetName((PetscObject)coordinates,"coordinates");
1187: VecGetArray(coordinates, &coords);
1188: for (k = 0; k < nVz; ++k) {
1189: PetscInt ind[3], d, off;
1191: ind[0] = 0;
1192: ind[1] = 0;
1193: ind[2] = k + da->zs;
1194: for (j = 0; j < nVy; ++j) {
1195: ind[1] = j + da->ys;
1196: for (i = 0; i < nVx; ++i) {
1197: const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;
1199: PetscSectionGetOffset(section, vertex, &off);
1200: ind[0] = i + da->xs;
1201: for (d = 0; d < dim; ++d) {
1202: coords[off+d] = h[d]*ind[d];
1203: }
1204: }
1205: }
1206: }
1207: VecRestoreArray(coordinates, &coords);
1208: DMSetCoordinateSection(dm, PETSC_DETERMINE, section);
1209: DMSetCoordinatesLocal(dm, coordinates);
1210: PetscSectionDestroy(§ion);
1211: VecDestroy(&coordinates);
1212: return(0);
1213: }
1217: PetscErrorCode DMProjectFunctionLocal_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1218: {
1219: PetscDS prob;
1220: PetscFE fe;
1221: PetscDualSpace sp;
1222: PetscQuadrature q;
1223: PetscSection section;
1224: PetscScalar *values;
1225: PetscInt numFields, numComp, numPoints, dim, dimEmbed, spDim, totDim, numValues, cStart, cEnd, f, c, v, d;
1226: PetscErrorCode ierr;
1229: DMGetDS(dm, &prob);
1230: PetscDSGetTotalDimension(prob, &totDim);
1231: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1232: PetscFEGetQuadrature(fe, &q);
1233: DMGetDefaultSection(dm, §ion);
1234: PetscSectionGetNumFields(section, &numFields);
1235: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1236: DMGetCoordinateDim(dm, &dimEmbed);
1237: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1238: DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
1239: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
1240: DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
1241: PetscQuadratureGetData(q, NULL, &numPoints, NULL, NULL);
1242: for (c = cStart; c < cEnd; ++c) {
1243: PetscFECellGeom geom;
1245: DMDAComputeCellGeometryFEM(dm, c, q, geom.v0, geom.J, NULL, &geom.detJ);
1246: geom.dim = dim;
1247: geom.dimEmbed = dimEmbed;
1248: for (f = 0, v = 0; f < numFields; ++f) {
1249: void * const ctx = ctxs ? ctxs[f] : NULL;
1251: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1252: PetscFEGetNumComponents(fe, &numComp);
1253: PetscFEGetDualSpace(fe, &sp);
1254: PetscDualSpaceGetDimension(sp, &spDim);
1255: for (d = 0; d < spDim; ++d) {
1256: PetscDualSpaceApply(sp, d, time, &geom, numComp, funcs[f], ctx, &values[v]);
1257: v += numComp;
1258: }
1259: }
1260: DMDAVecSetClosure(dm, section, localX, c, values, mode);
1261: }
1262: DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
1263: return(0);
1264: }
1268: PetscErrorCode DMComputeL2Diff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1269: {
1270: const PetscInt debug = 0;
1271: PetscDS prob;
1272: PetscFE fe;
1273: PetscQuadrature quad;
1274: PetscSection section;
1275: Vec localX;
1276: PetscScalar *funcVal;
1277: PetscReal *coords, *v0, *J, *invJ, detJ;
1278: PetscReal localDiff = 0.0;
1279: PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1280: PetscErrorCode ierr;
1283: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1284: DMGetDS(dm, &prob);
1285: PetscDSGetTotalComponents(prob, &Nc);
1286: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1287: PetscFEGetQuadrature(fe, &quad);
1288: DMGetDefaultSection(dm, §ion);
1289: PetscSectionGetNumFields(section, &numFields);
1290: DMGetLocalVector(dm, &localX);
1291: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1292: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1293: /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1294: PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1295: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1296: for (c = cStart; c < cEnd; ++c) {
1297: PetscScalar *x = NULL;
1298: PetscReal elemDiff = 0.0;
1300: DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);
1301: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1302: DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);
1304: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1305: void * const ctx = ctxs ? ctxs[field] : NULL;
1306: const PetscReal *quadPoints, *quadWeights;
1307: PetscReal *basis;
1308: PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;
1310: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1311: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1312: PetscFEGetDimension(fe, &numBasisFuncs);
1313: PetscFEGetNumComponents(fe, &numBasisComps);
1314: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
1315: if (debug) {
1316: char title[1024];
1317: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1318: DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
1319: }
1320: for (q = 0; q < numQuadPoints; ++q) {
1321: for (d = 0; d < dim; d++) {
1322: coords[d] = v0[d];
1323: for (e = 0; e < dim; e++) {
1324: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1325: }
1326: }
1327: (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);
1328: for (fc = 0; fc < numBasisComps; ++fc) {
1329: PetscScalar interpolant = 0.0;
1331: for (f = 0; f < numBasisFuncs; ++f) {
1332: const PetscInt fidx = f*numBasisComps+fc;
1333: interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
1334: }
1335: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
1336: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
1337: }
1338: }
1339: comp += numBasisComps;
1340: fieldOffset += numBasisFuncs*numBasisComps;
1341: }
1342: DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1343: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
1344: localDiff += elemDiff;
1345: }
1346: PetscFree5(funcVal,coords,v0,J,invJ);
1347: DMRestoreLocalVector(dm, &localX);
1348: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1349: *diff = PetscSqrtReal(*diff);
1350: return(0);
1351: }
1355: PetscErrorCode DMComputeL2GradientDiff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
1356: {
1357: const PetscInt debug = 0;
1358: PetscDS prob;
1359: PetscFE fe;
1360: PetscQuadrature quad;
1361: PetscSection section;
1362: Vec localX;
1363: PetscScalar *funcVal, *interpolantVec;
1364: PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
1365: PetscReal localDiff = 0.0;
1366: PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1367: PetscErrorCode ierr;
1370: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1371: DMGetDS(dm, &prob);
1372: PetscDSGetTotalComponents(prob, &Nc);
1373: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1374: PetscFEGetQuadrature(fe, &quad);
1375: DMGetDefaultSection(dm, §ion);
1376: PetscSectionGetNumFields(section, &numFields);
1377: DMGetLocalVector(dm, &localX);
1378: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1379: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1380: /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1381: PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
1382: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1383: for (c = cStart; c < cEnd; ++c) {
1384: PetscScalar *x = NULL;
1385: PetscReal elemDiff = 0.0;
1387: DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);
1388: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1389: DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);
1391: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1392: void * const ctx = ctxs ? ctxs[field] : NULL;
1393: const PetscReal *quadPoints, *quadWeights;
1394: PetscReal *basisDer;
1395: PetscInt Nq, Nb, Nc, q, d, e, fc, f, g;
1397: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1398: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
1399: PetscFEGetDimension(fe, &Nb);
1400: PetscFEGetNumComponents(fe, &Nc);
1401: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1402: if (debug) {
1403: char title[1024];
1404: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1405: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1406: }
1407: for (q = 0; q < Nq; ++q) {
1408: for (d = 0; d < dim; d++) {
1409: coords[d] = v0[d];
1410: for (e = 0; e < dim; e++) {
1411: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1412: }
1413: }
1414: (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);
1415: for (fc = 0; fc < Nc; ++fc) {
1416: PetscScalar interpolant = 0.0;
1418: for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
1419: for (f = 0; f < Nb; ++f) {
1420: const PetscInt fidx = f*Nc+fc;
1422: for (d = 0; d < dim; ++d) {
1423: realSpaceDer[d] = 0.0;
1424: for (g = 0; g < dim; ++g) {
1425: realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Nc+fidx)*dim+g];
1426: }
1427: interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
1428: }
1429: }
1430: for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
1431: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
1432: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
1433: }
1434: }
1435: comp += Nc;
1436: fieldOffset += Nb*Nc;
1437: }
1438: DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1439: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
1440: localDiff += elemDiff;
1441: }
1442: PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
1443: DMRestoreLocalVector(dm, &localX);
1444: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1445: *diff = PetscSqrtReal(*diff);
1446: return(0);
1447: }
1449: /* ------------------------------------------------------------------- */
1453: /*@C
1454: DMDAGetArray - Gets a work array for a DMDA
1456: Input Parameter:
1457: + da - information about my local patch
1458: - ghosted - do you want arrays for the ghosted or nonghosted patch
1460: Output Parameters:
1461: . vptr - array data structured
1463: Note: The vector values are NOT initialized and may have garbage in them, so you may need
1464: to zero them.
1466: Level: advanced
1468: .seealso: DMDARestoreArray()
1470: @*/
1471: PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
1472: {
1474: PetscInt j,i,xs,ys,xm,ym,zs,zm;
1475: char *iarray_start;
1476: void **iptr = (void**)vptr;
1477: DM_DA *dd = (DM_DA*)da->data;
1481: if (ghosted) {
1482: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1483: if (dd->arrayghostedin[i]) {
1484: *iptr = dd->arrayghostedin[i];
1485: iarray_start = (char*)dd->startghostedin[i];
1486: dd->arrayghostedin[i] = NULL;
1487: dd->startghostedin[i] = NULL;
1489: goto done;
1490: }
1491: }
1492: xs = dd->Xs;
1493: ys = dd->Ys;
1494: zs = dd->Zs;
1495: xm = dd->Xe-dd->Xs;
1496: ym = dd->Ye-dd->Ys;
1497: zm = dd->Ze-dd->Zs;
1498: } else {
1499: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1500: if (dd->arrayin[i]) {
1501: *iptr = dd->arrayin[i];
1502: iarray_start = (char*)dd->startin[i];
1503: dd->arrayin[i] = NULL;
1504: dd->startin[i] = NULL;
1506: goto done;
1507: }
1508: }
1509: xs = dd->xs;
1510: ys = dd->ys;
1511: zs = dd->zs;
1512: xm = dd->xe-dd->xs;
1513: ym = dd->ye-dd->ys;
1514: zm = dd->ze-dd->zs;
1515: }
1517: switch (da->dim) {
1518: case 1: {
1519: void *ptr;
1521: PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);
1523: ptr = (void*)(iarray_start - xs*sizeof(PetscScalar));
1524: *iptr = (void*)ptr;
1525: break;
1526: }
1527: case 2: {
1528: void **ptr;
1530: PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);
1532: ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
1533: for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
1534: *iptr = (void*)ptr;
1535: break;
1536: }
1537: case 3: {
1538: void ***ptr,**bptr;
1540: PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);
1542: ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1543: bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1544: for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1545: for (i=zs; i<zs+zm; i++) {
1546: for (j=ys; j<ys+ym; j++) {
1547: ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1548: }
1549: }
1550: *iptr = (void*)ptr;
1551: break;
1552: }
1553: default:
1554: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
1555: }
1557: done:
1558: /* add arrays to the checked out list */
1559: if (ghosted) {
1560: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1561: if (!dd->arrayghostedout[i]) {
1562: dd->arrayghostedout[i] = *iptr;
1563: dd->startghostedout[i] = iarray_start;
1564: break;
1565: }
1566: }
1567: } else {
1568: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1569: if (!dd->arrayout[i]) {
1570: dd->arrayout[i] = *iptr;
1571: dd->startout[i] = iarray_start;
1572: break;
1573: }
1574: }
1575: }
1576: return(0);
1577: }
1581: /*@C
1582: DMDARestoreArray - Restores an array of derivative types for a DMDA
1584: Input Parameter:
1585: + da - information about my local patch
1586: . ghosted - do you want arrays for the ghosted or nonghosted patch
1587: - vptr - array data structured to be passed to ad_FormFunctionLocal()
1589: Level: advanced
1591: .seealso: DMDAGetArray()
1593: @*/
1594: PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1595: {
1596: PetscInt i;
1597: void **iptr = (void**)vptr,*iarray_start = 0;
1598: DM_DA *dd = (DM_DA*)da->data;
1602: if (ghosted) {
1603: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1604: if (dd->arrayghostedout[i] == *iptr) {
1605: iarray_start = dd->startghostedout[i];
1606: dd->arrayghostedout[i] = NULL;
1607: dd->startghostedout[i] = NULL;
1608: break;
1609: }
1610: }
1611: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1612: if (!dd->arrayghostedin[i]) {
1613: dd->arrayghostedin[i] = *iptr;
1614: dd->startghostedin[i] = iarray_start;
1615: break;
1616: }
1617: }
1618: } else {
1619: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1620: if (dd->arrayout[i] == *iptr) {
1621: iarray_start = dd->startout[i];
1622: dd->arrayout[i] = NULL;
1623: dd->startout[i] = NULL;
1624: break;
1625: }
1626: }
1627: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1628: if (!dd->arrayin[i]) {
1629: dd->arrayin[i] = *iptr;
1630: dd->startin[i] = iarray_start;
1631: break;
1632: }
1633: }
1634: }
1635: return(0);
1636: }