Actual source code: dalocal.c
petsc-3.10.5 2019-03-28
2: /*
3: Code for manipulating distributed regular arrays in parallel.
4: */
6: #include <petsc/private/dmdaimpl.h>
7: #include <petscbt.h>
8: #include <petscsf.h>
9: #include <petscds.h>
10: #include <petscfe.h>
12: /*
13: This allows the DMDA vectors to properly tell MATLAB their dimensions
14: */
15: #if defined(PETSC_HAVE_MATLAB_ENGINE)
16: #include <engine.h> /* MATLAB include file */
17: #include <mex.h> /* MATLAB include file */
18: static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
19: {
21: PetscInt n,m;
22: Vec vec = (Vec)obj;
23: PetscScalar *array;
24: mxArray *mat;
25: DM da;
28: VecGetDM(vec, &da);
29: if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA");
30: DMDAGetGhostCorners(da,0,0,0,&m,&n,0);
32: VecGetArray(vec,&array);
33: #if !defined(PETSC_USE_COMPLEX)
34: mat = mxCreateDoubleMatrix(m,n,mxREAL);
35: #else
36: mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
37: #endif
38: PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));
39: PetscObjectName(obj);
40: engPutVariable((Engine*)mengine,obj->name,mat);
42: VecRestoreArray(vec,&array);
43: return(0);
44: }
45: #endif
48: PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g)
49: {
51: DM_DA *dd = (DM_DA*)da->data;
56: if (da->defaultSection) {
57: DMCreateLocalVector_Section_Private(da,g);
58: } else {
59: VecCreate(PETSC_COMM_SELF,g);
60: VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);
61: VecSetBlockSize(*g,dd->w);
62: VecSetType(*g,da->vectype);
63: VecSetDM(*g, da);
64: #if defined(PETSC_HAVE_MATLAB_ENGINE)
65: if (dd->w == 1 && da->dim == 2) {
66: PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);
67: }
68: #endif
69: }
70: return(0);
71: }
73: /*@
74: DMDAGetNumCells - Get the number of cells in the local piece of the DMDA. This includes ghost cells.
76: Input Parameter:
77: . dm - The DM object
79: Output Parameters:
80: + numCellsX - The number of local cells in the x-direction
81: . numCellsY - The number of local cells in the y-direction
82: . numCellsZ - The number of local cells in the z-direction
83: - numCells - The number of local cells
85: Level: developer
87: .seealso: DMDAGetCellPoint()
88: @*/
89: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
90: {
91: DM_DA *da = (DM_DA*) dm->data;
92: const PetscInt dim = dm->dim;
93: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
94: const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
98: if (numCellsX) {
100: *numCellsX = mx;
101: }
102: if (numCellsY) {
104: *numCellsY = my;
105: }
106: if (numCellsZ) {
108: *numCellsZ = mz;
109: }
110: if (numCells) {
112: *numCells = nC;
113: }
114: return(0);
115: }
117: /*@
118: DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the DMDA
120: Input Parameters:
121: + dm - The DM object
122: - i,j,k - The global indices for the cell
124: Output Parameters:
125: . point - The local DM point
127: Level: developer
129: .seealso: DMDAGetNumCells()
130: @*/
131: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
132: {
133: const PetscInt dim = dm->dim;
134: DMDALocalInfo info;
140: DMDAGetLocalInfo(dm, &info);
141: if (dim > 0) {if ((i < info.gxs) || (i >= info.gxs+info.gxm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "X index %d not in [%d, %d)", i, info.gxs, info.gxs+info.gxm);}
142: if (dim > 1) {if ((j < info.gys) || (j >= info.gys+info.gym)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %d not in [%d, %d)", j, info.gys, info.gys+info.gym);}
143: if (dim > 2) {if ((k < info.gzs) || (k >= info.gzs+info.gzm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %d not in [%d, %d)", k, info.gzs, info.gzs+info.gzm);}
144: *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0);
145: return(0);
146: }
148: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
149: {
150: DM_DA *da = (DM_DA*) dm->data;
151: const PetscInt dim = dm->dim;
152: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
153: const PetscInt nVx = mx+1;
154: const PetscInt nVy = dim > 1 ? (my+1) : 1;
155: const PetscInt nVz = dim > 2 ? (mz+1) : 1;
156: const PetscInt nV = nVx*nVy*nVz;
159: if (numVerticesX) {
161: *numVerticesX = nVx;
162: }
163: if (numVerticesY) {
165: *numVerticesY = nVy;
166: }
167: if (numVerticesZ) {
169: *numVerticesZ = nVz;
170: }
171: if (numVertices) {
173: *numVertices = nV;
174: }
175: return(0);
176: }
178: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
179: {
180: DM_DA *da = (DM_DA*) dm->data;
181: const PetscInt dim = dm->dim;
182: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
183: const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
184: const PetscInt nXF = (mx+1)*nxF;
185: const PetscInt nyF = mx*(dim > 2 ? mz : 1);
186: const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0;
187: const PetscInt nzF = mx*(dim > 1 ? my : 0);
188: const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0;
191: if (numXFacesX) {
193: *numXFacesX = nxF;
194: }
195: if (numXFaces) {
197: *numXFaces = nXF;
198: }
199: if (numYFacesY) {
201: *numYFacesY = nyF;
202: }
203: if (numYFaces) {
205: *numYFaces = nYF;
206: }
207: if (numZFacesZ) {
209: *numZFacesZ = nzF;
210: }
211: if (numZFaces) {
213: *numZFaces = nZF;
214: }
215: return(0);
216: }
218: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
219: {
220: const PetscInt dim = dm->dim;
221: PetscInt nC, nV, nXF, nYF, nZF;
227: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
228: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
229: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
230: if (height == 0) {
231: /* Cells */
232: if (pStart) *pStart = 0;
233: if (pEnd) *pEnd = nC;
234: } else if (height == 1) {
235: /* Faces */
236: if (pStart) *pStart = nC+nV;
237: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
238: } else if (height == dim) {
239: /* Vertices */
240: if (pStart) *pStart = nC;
241: if (pEnd) *pEnd = nC+nV;
242: } else if (height < 0) {
243: /* All points */
244: if (pStart) *pStart = 0;
245: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
246: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
247: return(0);
248: }
250: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
251: {
252: const PetscInt dim = dm->dim;
253: PetscInt nC, nV, nXF, nYF, nZF;
259: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
260: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
261: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
262: if (depth == dim) {
263: /* Cells */
264: if (pStart) *pStart = 0;
265: if (pEnd) *pEnd = nC;
266: } else if (depth == dim-1) {
267: /* Faces */
268: if (pStart) *pStart = nC+nV;
269: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
270: } else if (depth == 0) {
271: /* Vertices */
272: if (pStart) *pStart = nC;
273: if (pEnd) *pEnd = nC+nV;
274: } else if (depth < 0) {
275: /* All points */
276: if (pStart) *pStart = 0;
277: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
278: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
279: return(0);
280: }
282: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
283: {
284: const PetscInt dim = dm->dim;
285: PetscInt nC, nV, nXF, nYF, nZF;
289: *coneSize = 0;
290: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
291: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
292: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
293: switch (dim) {
294: case 2:
295: if (p >= 0) {
296: if (p < nC) {
297: *coneSize = 4;
298: } else if (p < nC+nV) {
299: *coneSize = 0;
300: } else if (p < nC+nV+nXF+nYF+nZF) {
301: *coneSize = 2;
302: } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
303: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
304: break;
305: case 3:
306: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
307: break;
308: }
309: return(0);
310: }
312: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
313: {
314: const PetscInt dim = dm->dim;
315: PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
319: if (!cone) {DMGetWorkArray(dm, 6, MPIU_INT, cone);}
320: DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);
321: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
322: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
323: switch (dim) {
324: case 2:
325: if (p >= 0) {
326: if (p < nC) {
327: const PetscInt cy = p / nCx;
328: const PetscInt cx = p % nCx;
330: (*cone)[0] = cy*nxF + cx + nC+nV;
331: (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
332: (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
333: (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
334: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
335: } else if (p < nC+nV) {
336: } else if (p < nC+nV+nXF) {
337: const PetscInt fy = (p - nC+nV) / nxF;
338: const PetscInt fx = (p - nC+nV) % nxF;
340: (*cone)[0] = fy*nVx + fx + nC;
341: (*cone)[1] = fy*nVx + fx + 1 + nC;
342: } else if (p < nC+nV+nXF+nYF) {
343: const PetscInt fx = (p - nC+nV+nXF) / nyF;
344: const PetscInt fy = (p - nC+nV+nXF) % nyF;
346: (*cone)[0] = fy*nVx + fx + nC;
347: (*cone)[1] = fy*nVx + fx + nVx + nC;
348: } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
349: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
350: break;
351: case 3:
352: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
353: break;
354: }
355: return(0);
356: }
358: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
359: {
363: DMGetWorkArray(dm, 6, MPIU_INT, cone);
364: return(0);
365: }
367: /*@C
368: DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with
369: different numbers of dofs on vertices, cells, and faces in each direction.
371: Input Parameters:
372: + dm- The DMDA
373: . numFields - The number of fields
374: . numComp - The number of components in each field
375: . numDof - The number of dofs per dimension for each field
376: . numFaceDof - The number of dofs per face for each field and direction, or NULL
378: Level: developer
380: Note:
381: The default DMDA numbering is as follows:
383: - Cells: [0, nC)
384: - Vertices: [nC, nC+nV)
385: - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir
386: - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir
387: - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir
389: We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment.
390: @*/
391: PetscErrorCode DMDACreateSection(DM dm, const PetscInt numComp[], const PetscInt numDof[], const PetscInt numFaceDof[], PetscSection *s)
392: {
393: PetscSection section;
394: const PetscInt dim = dm->dim;
395: PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0};
396: PetscBT isLeaf;
397: PetscSF sf;
398: PetscMPIInt rank;
399: const PetscMPIInt *neighbors;
400: PetscInt *localPoints;
401: PetscSFNode *remotePoints;
402: PetscInt nleaves = 0, nleavesCheck = 0, nL = 0;
403: PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF;
404: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd;
405: PetscInt f, v, c, xf, yf, zf, xn, yn, zn;
406: PetscErrorCode ierr;
411: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
412: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
413: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
414: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
415: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
416: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
417: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
418: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
419: xfStart = vEnd; xfEnd = xfStart+nXF;
420: yfStart = xfEnd; yfEnd = yfStart+nYF;
421: zfStart = yfEnd; zfEnd = zfStart+nZF;
422: if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd);
423: /* Create local section */
424: DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);
425: for (f = 0; f < numFields; ++f) {
426: numVertexTotDof += numDof[f*(dim+1)+0];
427: numCellTotDof += numDof[f*(dim+1)+dim];
428: numFaceTotDof[0] += dim > 0 ? (numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]) : 0;
429: numFaceTotDof[1] += dim > 1 ? (numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]) : 0;
430: numFaceTotDof[2] += dim > 2 ? (numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]) : 0;
431: }
432: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
433: if (numFields > 0) {
434: PetscSectionSetNumFields(section, numFields);
435: for (f = 0; f < numFields; ++f) {
436: const char *name;
438: DMDAGetFieldName(dm, f, &name);
439: PetscSectionSetFieldName(section, f, name ? name : "Field");
440: if (numComp) {
441: PetscSectionSetFieldComponents(section, f, numComp[f]);
442: }
443: }
444: }
445: PetscSectionSetChart(section, pStart, pEnd);
446: for (v = vStart; v < vEnd; ++v) {
447: for (f = 0; f < numFields; ++f) {
448: PetscSectionSetFieldDof(section, v, f, numDof[f*(dim+1)+0]);
449: }
450: PetscSectionSetDof(section, v, numVertexTotDof);
451: }
452: for (xf = xfStart; xf < xfEnd; ++xf) {
453: for (f = 0; f < numFields; ++f) {
454: PetscSectionSetFieldDof(section, xf, f, numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]);
455: }
456: PetscSectionSetDof(section, xf, numFaceTotDof[0]);
457: }
458: for (yf = yfStart; yf < yfEnd; ++yf) {
459: for (f = 0; f < numFields; ++f) {
460: PetscSectionSetFieldDof(section, yf, f, numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]);
461: }
462: PetscSectionSetDof(section, yf, numFaceTotDof[1]);
463: }
464: for (zf = zfStart; zf < zfEnd; ++zf) {
465: for (f = 0; f < numFields; ++f) {
466: PetscSectionSetFieldDof(section, zf, f, numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]);
467: }
468: PetscSectionSetDof(section, zf, numFaceTotDof[2]);
469: }
470: for (c = cStart; c < cEnd; ++c) {
471: for (f = 0; f < numFields; ++f) {
472: PetscSectionSetFieldDof(section, c, f, numDof[f*(dim+1)+dim]);
473: }
474: PetscSectionSetDof(section, c, numCellTotDof);
475: }
476: PetscSectionSetUp(section);
477: /* Create mesh point SF */
478: PetscBTCreate(pEnd-pStart, &isLeaf);
479: DMDAGetNeighbors(dm, &neighbors);
480: for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
481: for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
482: for (xn = 0; xn < 3; ++xn) {
483: const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
484: const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
485: PetscInt xv, yv, zv;
487: if (neighbor >= 0 && neighbor < rank) {
488: if (xp < 0) { /* left */
489: if (yp < 0) { /* bottom */
490: if (zp < 0) { /* back */
491: const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* left bottom back vertex */
492: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
493: } else if (zp > 0) { /* front */
494: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* left bottom front vertex */
495: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
496: } else {
497: for (zv = 0; zv < nVz; ++zv) {
498: const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; /* left bottom vertices */
499: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
500: }
501: }
502: } else if (yp > 0) { /* top */
503: if (zp < 0) { /* back */
504: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* left top back vertex */
505: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
506: } else if (zp > 0) { /* front */
507: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* left top front vertex */
508: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
509: } else {
510: for (zv = 0; zv < nVz; ++zv) {
511: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* left top vertices */
512: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
513: }
514: }
515: } else {
516: if (zp < 0) { /* back */
517: for (yv = 0; yv < nVy; ++yv) {
518: const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* left back vertices */
519: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
520: }
521: } else if (zp > 0) { /* front */
522: for (yv = 0; yv < nVy; ++yv) {
523: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* left front vertices */
524: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
525: }
526: } else {
527: for (zv = 0; zv < nVz; ++zv) {
528: for (yv = 0; yv < nVy; ++yv) {
529: const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; /* left vertices */
530: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
531: }
532: }
533: #if 0
534: for (xf = 0; xf < nxF; ++xf) {
535: /* THIS IS WRONG */
536: const PetscInt localFace = 0 + nC+nV; /* left faces */
537: if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
538: }
539: #endif
540: }
541: }
542: } else if (xp > 0) { /* right */
543: if (yp < 0) { /* bottom */
544: if (zp < 0) { /* back */
545: const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* right bottom back vertex */
546: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
547: } else if (zp > 0) { /* front */
548: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* right bottom front vertex */
549: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
550: } else {
551: for (zv = 0; zv < nVz; ++zv) {
552: const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* right bottom vertices */
553: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
554: }
555: }
556: } else if (yp > 0) { /* top */
557: if (zp < 0) { /* back */
558: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
559: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
560: } else if (zp > 0) { /* front */
561: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
562: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
563: } else {
564: for (zv = 0; zv < nVz; ++zv) {
565: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
566: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
567: }
568: }
569: } else {
570: if (zp < 0) { /* back */
571: for (yv = 0; yv < nVy; ++yv) {
572: const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
573: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
574: }
575: } else if (zp > 0) { /* front */
576: for (yv = 0; yv < nVy; ++yv) {
577: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
578: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
579: }
580: } else {
581: for (zv = 0; zv < nVz; ++zv) {
582: for (yv = 0; yv < nVy; ++yv) {
583: const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
584: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
585: }
586: }
587: #if 0
588: for (xf = 0; xf < nxF; ++xf) {
589: /* THIS IS WRONG */
590: const PetscInt localFace = 0 + nC+nV; /* right faces */
591: if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
592: }
593: #endif
594: }
595: }
596: } else {
597: if (yp < 0) { /* bottom */
598: if (zp < 0) { /* back */
599: for (xv = 0; xv < nVx; ++xv) {
600: const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; /* bottom back vertices */
601: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
602: }
603: } else if (zp > 0) { /* front */
604: for (xv = 0; xv < nVx; ++xv) {
605: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* bottom front vertices */
606: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
607: }
608: } else {
609: for (zv = 0; zv < nVz; ++zv) {
610: for (xv = 0; xv < nVx; ++xv) {
611: const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; /* bottom vertices */
612: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
613: }
614: }
615: #if 0
616: for (yf = 0; yf < nyF; ++yf) {
617: /* THIS IS WRONG */
618: const PetscInt localFace = 0 + nC+nV; /* bottom faces */
619: if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves;
620: }
621: #endif
622: }
623: } else if (yp > 0) { /* top */
624: if (zp < 0) { /* back */
625: for (xv = 0; xv < nVx; ++xv) {
626: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
627: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
628: }
629: } else if (zp > 0) { /* front */
630: for (xv = 0; xv < nVx; ++xv) {
631: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
632: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
633: }
634: } else {
635: for (zv = 0; zv < nVz; ++zv) {
636: for (xv = 0; xv < nVx; ++xv) {
637: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
638: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
639: }
640: }
641: #if 0
642: for (yf = 0; yf < nyF; ++yf) {
643: /* THIS IS WRONG */
644: const PetscInt localFace = 0 + nC+nV; /* top faces */
645: if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves;
646: }
647: #endif
648: }
649: } else {
650: if (zp < 0) { /* back */
651: for (yv = 0; yv < nVy; ++yv) {
652: for (xv = 0; xv < nVx; ++xv) {
653: const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; /* back vertices */
654: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
655: }
656: }
657: #if 0
658: for (zf = 0; zf < nzF; ++zf) {
659: /* THIS IS WRONG */
660: const PetscInt localFace = 0 + nC+nV; /* back faces */
661: if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
662: }
663: #endif
664: } else if (zp > 0) { /* front */
665: for (yv = 0; yv < nVy; ++yv) {
666: for (xv = 0; xv < nVx; ++xv) {
667: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
668: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
669: }
670: }
671: #if 0
672: for (zf = 0; zf < nzF; ++zf) {
673: /* THIS IS WRONG */
674: const PetscInt localFace = 0 + nC+nV; /* front faces */
675: if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
676: }
677: #endif
678: } else {
679: /* Nothing is shared from the interior */
680: }
681: }
682: }
683: }
684: }
685: }
686: }
687: PetscBTMemzero(pEnd-pStart, isLeaf);
688: PetscMalloc1(nleaves,&localPoints);
689: PetscMalloc1(nleaves,&remotePoints);
690: for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
691: for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
692: for (xn = 0; xn < 3; ++xn) {
693: const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
694: const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
695: PetscInt xv, yv, zv;
697: if (neighbor >= 0 && neighbor < rank) {
698: if (xp < 0) { /* left */
699: if (yp < 0) { /* bottom */
700: if (zp < 0) { /* back */
701: const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* left bottom back vertex */
702: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
704: if (!PetscBTLookupSet(isLeaf, localVertex)) {
705: localPoints[nL] = localVertex;
706: remotePoints[nL].rank = neighbor;
707: remotePoints[nL].index = remoteVertex;
708: ++nL;
709: }
710: } else if (zp > 0) { /* front */
711: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* left bottom front vertex */
712: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
714: if (!PetscBTLookupSet(isLeaf, localVertex)) {
715: localPoints[nL] = localVertex;
716: remotePoints[nL].rank = neighbor;
717: remotePoints[nL].index = remoteVertex;
718: ++nL;
719: }
720: } else {
721: for (zv = 0; zv < nVz; ++zv) {
722: const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; /* left bottom vertices */
723: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
725: if (!PetscBTLookupSet(isLeaf, localVertex)) {
726: localPoints[nL] = localVertex;
727: remotePoints[nL].rank = neighbor;
728: remotePoints[nL].index = remoteVertex;
729: ++nL;
730: }
731: }
732: }
733: } else if (yp > 0) { /* top */
734: if (zp < 0) { /* back */
735: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* left top back vertex */
736: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
738: if (!PetscBTLookupSet(isLeaf, localVertex)) {
739: localPoints[nL] = localVertex;
740: remotePoints[nL].rank = neighbor;
741: remotePoints[nL].index = remoteVertex;
742: ++nL;
743: }
744: } else if (zp > 0) { /* front */
745: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* left top front vertex */
746: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
748: if (!PetscBTLookupSet(isLeaf, localVertex)) {
749: localPoints[nL] = localVertex;
750: remotePoints[nL].rank = neighbor;
751: remotePoints[nL].index = remoteVertex;
752: ++nL;
753: }
754: } else {
755: for (zv = 0; zv < nVz; ++zv) {
756: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* left top vertices */
757: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
759: if (!PetscBTLookupSet(isLeaf, localVertex)) {
760: localPoints[nL] = localVertex;
761: remotePoints[nL].rank = neighbor;
762: remotePoints[nL].index = remoteVertex;
763: ++nL;
764: }
765: }
766: }
767: } else {
768: if (zp < 0) { /* back */
769: for (yv = 0; yv < nVy; ++yv) {
770: const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* left back vertices */
771: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
773: if (!PetscBTLookupSet(isLeaf, localVertex)) {
774: localPoints[nL] = localVertex;
775: remotePoints[nL].rank = neighbor;
776: remotePoints[nL].index = remoteVertex;
777: ++nL;
778: }
779: }
780: } else if (zp > 0) { /* front */
781: for (yv = 0; yv < nVy; ++yv) {
782: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* left front vertices */
783: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
785: if (!PetscBTLookupSet(isLeaf, localVertex)) {
786: localPoints[nL] = localVertex;
787: remotePoints[nL].rank = neighbor;
788: remotePoints[nL].index = remoteVertex;
789: ++nL;
790: }
791: }
792: } else {
793: for (zv = 0; zv < nVz; ++zv) {
794: for (yv = 0; yv < nVy; ++yv) {
795: const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; /* left vertices */
796: const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
798: if (!PetscBTLookupSet(isLeaf, localVertex)) {
799: localPoints[nL] = localVertex;
800: remotePoints[nL].rank = neighbor;
801: remotePoints[nL].index = remoteVertex;
802: ++nL;
803: }
804: }
805: }
806: #if 0
807: for (xf = 0; xf < nxF; ++xf) {
808: /* THIS IS WRONG */
809: const PetscInt localFace = 0 + nC+nV; /* left faces */
810: const PetscInt remoteFace = 0 + nC+nV;
812: if (!PetscBTLookupSet(isLeaf, localFace)) {
813: localPoints[nL] = localFace;
814: remotePoints[nL].rank = neighbor;
815: remotePoints[nL].index = remoteFace;
816: }
817: }
818: #endif
819: }
820: }
821: } else if (xp > 0) { /* right */
822: if (yp < 0) { /* bottom */
823: if (zp < 0) { /* back */
824: const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* right bottom back vertex */
825: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
827: if (!PetscBTLookupSet(isLeaf, localVertex)) {
828: localPoints[nL] = localVertex;
829: remotePoints[nL].rank = neighbor;
830: remotePoints[nL].index = remoteVertex;
831: ++nL;
832: }
833: } else if (zp > 0) { /* front */
834: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* right bottom front vertex */
835: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
837: if (!PetscBTLookupSet(isLeaf, localVertex)) {
838: localPoints[nL] = localVertex;
839: remotePoints[nL].rank = neighbor;
840: remotePoints[nL].index = remoteVertex;
841: ++nL;
842: }
843: } else {
844: nleavesCheck += nVz;
845: for (zv = 0; zv < nVz; ++zv) {
846: const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* right bottom vertices */
847: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
849: if (!PetscBTLookupSet(isLeaf, localVertex)) {
850: localPoints[nL] = localVertex;
851: remotePoints[nL].rank = neighbor;
852: remotePoints[nL].index = remoteVertex;
853: ++nL;
854: }
855: }
856: }
857: } else if (yp > 0) { /* top */
858: if (zp < 0) { /* back */
859: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
860: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
862: if (!PetscBTLookupSet(isLeaf, localVertex)) {
863: localPoints[nL] = localVertex;
864: remotePoints[nL].rank = neighbor;
865: remotePoints[nL].index = remoteVertex;
866: ++nL;
867: }
868: } else if (zp > 0) { /* front */
869: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
870: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
872: if (!PetscBTLookupSet(isLeaf, localVertex)) {
873: localPoints[nL] = localVertex;
874: remotePoints[nL].rank = neighbor;
875: remotePoints[nL].index = remoteVertex;
876: ++nL;
877: }
878: } else {
879: for (zv = 0; zv < nVz; ++zv) {
880: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
881: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
883: if (!PetscBTLookupSet(isLeaf, localVertex)) {
884: localPoints[nL] = localVertex;
885: remotePoints[nL].rank = neighbor;
886: remotePoints[nL].index = remoteVertex;
887: ++nL;
888: }
889: }
890: }
891: } else {
892: if (zp < 0) { /* back */
893: for (yv = 0; yv < nVy; ++yv) {
894: const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
895: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
897: if (!PetscBTLookupSet(isLeaf, localVertex)) {
898: localPoints[nL] = localVertex;
899: remotePoints[nL].rank = neighbor;
900: remotePoints[nL].index = remoteVertex;
901: ++nL;
902: }
903: }
904: } else if (zp > 0) { /* front */
905: for (yv = 0; yv < nVy; ++yv) {
906: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
907: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
909: if (!PetscBTLookupSet(isLeaf, localVertex)) {
910: localPoints[nL] = localVertex;
911: remotePoints[nL].rank = neighbor;
912: remotePoints[nL].index = remoteVertex;
913: ++nL;
914: }
915: }
916: } else {
917: for (zv = 0; zv < nVz; ++zv) {
918: for (yv = 0; yv < nVy; ++yv) {
919: const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
920: const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
922: if (!PetscBTLookupSet(isLeaf, localVertex)) {
923: localPoints[nL] = localVertex;
924: remotePoints[nL].rank = neighbor;
925: remotePoints[nL].index = remoteVertex;
926: ++nL;
927: }
928: }
929: }
930: #if 0
931: for (xf = 0; xf < nxF; ++xf) {
932: /* THIS IS WRONG */
933: const PetscInt localFace = 0 + nC+nV; /* right faces */
934: const PetscInt remoteFace = 0 + nC+nV;
936: if (!PetscBTLookupSet(isLeaf, localFace)) {
937: localPoints[nL] = localFace;
938: remotePoints[nL].rank = neighbor;
939: remotePoints[nL].index = remoteFace;
940: ++nL;
941: }
942: }
943: #endif
944: }
945: }
946: } else {
947: if (yp < 0) { /* bottom */
948: if (zp < 0) { /* back */
949: for (xv = 0; xv < nVx; ++xv) {
950: const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; /* bottom back vertices */
951: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
953: if (!PetscBTLookupSet(isLeaf, localVertex)) {
954: localPoints[nL] = localVertex;
955: remotePoints[nL].rank = neighbor;
956: remotePoints[nL].index = remoteVertex;
957: ++nL;
958: }
959: }
960: } else if (zp > 0) { /* front */
961: for (xv = 0; xv < nVx; ++xv) {
962: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* bottom front vertices */
963: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
965: if (!PetscBTLookupSet(isLeaf, localVertex)) {
966: localPoints[nL] = localVertex;
967: remotePoints[nL].rank = neighbor;
968: remotePoints[nL].index = remoteVertex;
969: ++nL;
970: }
971: }
972: } else {
973: for (zv = 0; zv < nVz; ++zv) {
974: for (xv = 0; xv < nVx; ++xv) {
975: const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; /* bottom vertices */
976: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
978: if (!PetscBTLookupSet(isLeaf, localVertex)) {
979: localPoints[nL] = localVertex;
980: remotePoints[nL].rank = neighbor;
981: remotePoints[nL].index = remoteVertex;
982: ++nL;
983: }
984: }
985: }
986: #if 0
987: for (yf = 0; yf < nyF; ++yf) {
988: /* THIS IS WRONG */
989: const PetscInt localFace = 0 + nC+nV; /* bottom faces */
990: const PetscInt remoteFace = 0 + nC+nV;
992: if (!PetscBTLookupSet(isLeaf, localFace)) {
993: localPoints[nL] = localFace;
994: remotePoints[nL].rank = neighbor;
995: remotePoints[nL].index = remoteFace;
996: ++nL;
997: }
998: }
999: #endif
1000: }
1001: } else if (yp > 0) { /* top */
1002: if (zp < 0) { /* back */
1003: for (xv = 0; xv < nVx; ++xv) {
1004: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
1005: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1007: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1008: localPoints[nL] = localVertex;
1009: remotePoints[nL].rank = neighbor;
1010: remotePoints[nL].index = remoteVertex;
1011: ++nL;
1012: }
1013: }
1014: } else if (zp > 0) { /* front */
1015: for (xv = 0; xv < nVx; ++xv) {
1016: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
1017: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1019: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1020: localPoints[nL] = localVertex;
1021: remotePoints[nL].rank = neighbor;
1022: remotePoints[nL].index = remoteVertex;
1023: ++nL;
1024: }
1025: }
1026: } else {
1027: for (zv = 0; zv < nVz; ++zv) {
1028: for (xv = 0; xv < nVx; ++xv) {
1029: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
1030: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1032: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1033: localPoints[nL] = localVertex;
1034: remotePoints[nL].rank = neighbor;
1035: remotePoints[nL].index = remoteVertex;
1036: ++nL;
1037: }
1038: }
1039: }
1040: #if 0
1041: for (yf = 0; yf < nyF; ++yf) {
1042: /* THIS IS WRONG */
1043: const PetscInt localFace = 0 + nC+nV; /* top faces */
1044: const PetscInt remoteFace = 0 + nC+nV;
1046: if (!PetscBTLookupSet(isLeaf, localFace)) {
1047: localPoints[nL] = localFace;
1048: remotePoints[nL].rank = neighbor;
1049: remotePoints[nL].index = remoteFace;
1050: ++nL;
1051: }
1052: }
1053: #endif
1054: }
1055: } else {
1056: if (zp < 0) { /* back */
1057: for (yv = 0; yv < nVy; ++yv) {
1058: for (xv = 0; xv < nVx; ++xv) {
1059: const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; /* back vertices */
1060: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1062: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1063: localPoints[nL] = localVertex;
1064: remotePoints[nL].rank = neighbor;
1065: remotePoints[nL].index = remoteVertex;
1066: ++nL;
1067: }
1068: }
1069: }
1070: #if 0
1071: for (zf = 0; zf < nzF; ++zf) {
1072: /* THIS IS WRONG */
1073: const PetscInt localFace = 0 + nC+nV; /* back faces */
1074: const PetscInt remoteFace = 0 + nC+nV;
1076: if (!PetscBTLookupSet(isLeaf, localFace)) {
1077: localPoints[nL] = localFace;
1078: remotePoints[nL].rank = neighbor;
1079: remotePoints[nL].index = remoteFace;
1080: ++nL;
1081: }
1082: }
1083: #endif
1084: } else if (zp > 0) { /* front */
1085: for (yv = 0; yv < nVy; ++yv) {
1086: for (xv = 0; xv < nVx; ++xv) {
1087: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
1088: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1090: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1091: localPoints[nL] = localVertex;
1092: remotePoints[nL].rank = neighbor;
1093: remotePoints[nL].index = remoteVertex;
1094: ++nL;
1095: }
1096: }
1097: }
1098: #if 0
1099: for (zf = 0; zf < nzF; ++zf) {
1100: /* THIS IS WRONG */
1101: const PetscInt localFace = 0 + nC+nV; /* front faces */
1102: const PetscInt remoteFace = 0 + nC+nV;
1104: if (!PetscBTLookupSet(isLeaf, localFace)) {
1105: localPoints[nL] = localFace;
1106: remotePoints[nL].rank = neighbor;
1107: remotePoints[nL].index = remoteFace;
1108: ++nL;
1109: }
1110: }
1111: #endif
1112: } else {
1113: /* Nothing is shared from the interior */
1114: }
1115: }
1116: }
1117: }
1118: }
1119: }
1120: }
1121: PetscBTDestroy(&isLeaf);
1122: /* Remove duplication in leaf determination */
1123: if (nleaves != nL) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck);
1124: PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);
1125: PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1126: DMSetPointSF(dm, sf);
1127: PetscSFDestroy(&sf);
1128: *s = section;
1129: return(0);
1130: }
1132: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
1133: {
1134: DM_DA *da = (DM_DA *) dm->data;
1135: Vec coordinates;
1136: PetscSection section;
1137: PetscScalar *coords;
1138: PetscReal h[3];
1139: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
1144: DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);
1145: if (dim > 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"The following code only works for dim <= 3");
1146: h[0] = (xu - xl)/M;
1147: h[1] = (yu - yl)/N;
1148: h[2] = (zu - zl)/P;
1149: DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
1150: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
1151: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);
1152: PetscSectionSetNumFields(section, 1);
1153: PetscSectionSetFieldComponents(section, 0, dim);
1154: PetscSectionSetChart(section, vStart, vEnd);
1155: for (v = vStart; v < vEnd; ++v) {
1156: PetscSectionSetDof(section, v, dim);
1157: }
1158: PetscSectionSetUp(section);
1159: PetscSectionGetStorageSize(section, &size);
1160: VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
1161: PetscObjectSetName((PetscObject)coordinates,"coordinates");
1162: VecGetArray(coordinates, &coords);
1163: for (k = 0; k < nVz; ++k) {
1164: PetscInt ind[3], d, off;
1166: ind[0] = 0;
1167: ind[1] = 0;
1168: ind[2] = k + da->zs;
1169: for (j = 0; j < nVy; ++j) {
1170: ind[1] = j + da->ys;
1171: for (i = 0; i < nVx; ++i) {
1172: const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;
1174: PetscSectionGetOffset(section, vertex, &off);
1175: ind[0] = i + da->xs;
1176: for (d = 0; d < dim; ++d) {
1177: coords[off+d] = h[d]*ind[d];
1178: }
1179: }
1180: }
1181: }
1182: VecRestoreArray(coordinates, &coords);
1183: DMSetCoordinateSection(dm, PETSC_DETERMINE, section);
1184: DMSetCoordinatesLocal(dm, coordinates);
1185: PetscSectionDestroy(§ion);
1186: VecDestroy(&coordinates);
1187: return(0);
1188: }
1190: PetscErrorCode DMProjectFunctionLocal_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1191: {
1192: PetscDS prob;
1193: PetscFE fe;
1194: PetscDualSpace sp;
1195: PetscQuadrature q;
1196: PetscSection section;
1197: PetscScalar *values;
1198: PetscInt numFields, Nc, dim, dimEmbed, spDim, totDim, numValues, cStart, cEnd, f, c, v, d;
1199: PetscFEGeom *geom;
1200: PetscErrorCode ierr;
1203: DMGetDS(dm, &prob);
1204: PetscDSGetTotalDimension(prob, &totDim);
1205: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1206: PetscFEGetQuadrature(fe, &q);
1207: DMGetSection(dm, §ion);
1208: PetscSectionGetNumFields(section, &numFields);
1209: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1210: DMGetCoordinateDim(dm, &dimEmbed);
1211: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1212: DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
1213: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
1214: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
1215: PetscFEGeomCreate(q,1,dimEmbed,PETSC_FALSE,&geom);
1216: for (c = cStart; c < cEnd; ++c) {
1217: DMDAComputeCellGeometryFEM(dm, c, q, geom->v, geom->J, NULL, geom->detJ);
1218: for (f = 0, v = 0; f < numFields; ++f) {
1219: void * const ctx = ctxs ? ctxs[f] : NULL;
1221: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1222: PetscFEGetNumComponents(fe, &Nc);
1223: PetscFEGetDualSpace(fe, &sp);
1224: PetscDualSpaceGetDimension(sp, &spDim);
1225: for (d = 0; d < spDim; ++d) {
1226: PetscDualSpaceApply(sp, d, time, geom, Nc, funcs[f], ctx, &values[v]);
1227: }
1228: }
1229: DMDAVecSetClosure(dm, section, localX, c, values, mode);
1230: }
1231: PetscFEGeomDestroy(&geom);
1232: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
1233: return(0);
1234: }
1236: PetscErrorCode DMComputeL2Diff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1237: {
1238: const PetscInt debug = 0;
1239: PetscDS prob;
1240: PetscFE fe;
1241: PetscQuadrature quad;
1242: PetscSection section;
1243: Vec localX;
1244: PetscScalar *funcVal;
1245: PetscReal *coords, *v0, *J, *invJ, detJ;
1246: PetscReal localDiff = 0.0;
1247: PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1248: PetscErrorCode ierr;
1251: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1252: DMGetDS(dm, &prob);
1253: PetscDSGetTotalComponents(prob, &Nc);
1254: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1255: PetscFEGetQuadrature(fe, &quad);
1256: DMGetSection(dm, §ion);
1257: PetscSectionGetNumFields(section, &numFields);
1258: DMGetLocalVector(dm, &localX);
1259: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1260: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1261: /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1262: PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1263: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1264: for (c = cStart; c < cEnd; ++c) {
1265: PetscScalar *x = NULL;
1266: PetscReal elemDiff = 0.0;
1268: DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);
1269: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1270: DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);
1272: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1273: void * const ctx = ctxs ? ctxs[field] : NULL;
1274: const PetscReal *quadPoints, *quadWeights;
1275: PetscReal *basis;
1276: PetscInt qNc, Nq, Nb, Nc, q, d, e, fc, f;
1278: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1279: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1280: PetscFEGetDimension(fe, &Nb);
1281: PetscFEGetNumComponents(fe, &Nc);
1282: if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Quadrature components %D != %D field components\n", qNc, Nc);
1283: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
1284: if (debug) {
1285: char title[1024];
1286: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1287: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1288: }
1289: for (q = 0; q < Nq; ++q) {
1290: for (d = 0; d < dim; d++) {
1291: coords[d] = v0[d];
1292: for (e = 0; e < dim; e++) {
1293: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1294: }
1295: }
1296: (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);
1297: for (fc = 0; fc < Nc; ++fc) {
1298: PetscScalar interpolant = 0.0;
1300: for (f = 0; f < Nb; ++f) {
1301: interpolant += x[fieldOffset+f]*basis[(q*Nb+f)*Nc+fc];
1302: }
1303: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+fc]*detJ);}
1304: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ;
1305: }
1306: }
1307: comp += Nc;
1308: fieldOffset += Nb;
1309: }
1310: DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1311: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
1312: localDiff += elemDiff;
1313: }
1314: PetscFree5(funcVal,coords,v0,J,invJ);
1315: DMRestoreLocalVector(dm, &localX);
1316: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1317: *diff = PetscSqrtReal(*diff);
1318: return(0);
1319: }
1321: PetscErrorCode DMComputeL2GradientDiff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
1322: {
1323: const PetscInt debug = 0;
1324: PetscDS prob;
1325: PetscFE fe;
1326: PetscQuadrature quad;
1327: PetscSection section;
1328: Vec localX;
1329: PetscScalar *funcVal, *interpolantVec;
1330: PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
1331: PetscReal localDiff = 0.0;
1332: PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1333: PetscErrorCode ierr;
1336: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1337: DMGetDS(dm, &prob);
1338: PetscDSGetTotalComponents(prob, &Nc);
1339: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1340: PetscFEGetQuadrature(fe, &quad);
1341: DMGetSection(dm, §ion);
1342: PetscSectionGetNumFields(section, &numFields);
1343: DMGetLocalVector(dm, &localX);
1344: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1345: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1346: /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1347: PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
1348: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1349: for (c = cStart; c < cEnd; ++c) {
1350: PetscScalar *x = NULL;
1351: PetscReal elemDiff = 0.0;
1353: DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);
1354: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1355: DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);
1357: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1358: void * const ctx = ctxs ? ctxs[field] : NULL;
1359: const PetscReal *quadPoints, *quadWeights;
1360: PetscReal *basisDer;
1361: PetscInt qNc, Nq, Nb, Nc, q, d, e, fc, f, g;
1363: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1364: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1365: PetscFEGetDimension(fe, &Nb);
1366: PetscFEGetNumComponents(fe, &Nc);
1367: if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Quadrature components %D != %D field components\n", qNc, Nc);
1368: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1369: if (debug) {
1370: char title[1024];
1371: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1372: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1373: }
1374: for (q = 0; q < Nq; ++q) {
1375: for (d = 0; d < dim; d++) {
1376: coords[d] = v0[d];
1377: for (e = 0; e < dim; e++) {
1378: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1379: }
1380: }
1381: (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);
1382: for (fc = 0; fc < Nc; ++fc) {
1383: PetscScalar interpolant = 0.0;
1385: for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
1386: for (f = 0; f < Nb; ++f) {
1387: for (d = 0; d < dim; ++d) {
1388: realSpaceDer[d] = 0.0;
1389: for (g = 0; g < dim; ++g) {
1390: realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb+f)*dim+g];
1391: }
1392: interpolantVec[d] += x[fieldOffset+f]*realSpaceDer[d];
1393: }
1394: }
1395: for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
1396: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ);}
1397: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ;
1398: }
1399: }
1400: comp += Nc;
1401: fieldOffset += Nb*Nc;
1402: }
1403: DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1404: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
1405: localDiff += elemDiff;
1406: }
1407: PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
1408: DMRestoreLocalVector(dm, &localX);
1409: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1410: *diff = PetscSqrtReal(*diff);
1411: return(0);
1412: }
1414: /* ------------------------------------------------------------------- */
1416: /*@C
1417: DMDAGetArray - Gets a work array for a DMDA
1419: Input Parameter:
1420: + da - information about my local patch
1421: - ghosted - do you want arrays for the ghosted or nonghosted patch
1423: Output Parameters:
1424: . vptr - array data structured
1426: Note: The vector values are NOT initialized and may have garbage in them, so you may need
1427: to zero them.
1429: Level: advanced
1431: .seealso: DMDARestoreArray()
1433: @*/
1434: PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
1435: {
1437: PetscInt j,i,xs,ys,xm,ym,zs,zm;
1438: char *iarray_start;
1439: void **iptr = (void**)vptr;
1440: DM_DA *dd = (DM_DA*)da->data;
1444: if (ghosted) {
1445: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1446: if (dd->arrayghostedin[i]) {
1447: *iptr = dd->arrayghostedin[i];
1448: iarray_start = (char*)dd->startghostedin[i];
1449: dd->arrayghostedin[i] = NULL;
1450: dd->startghostedin[i] = NULL;
1452: goto done;
1453: }
1454: }
1455: xs = dd->Xs;
1456: ys = dd->Ys;
1457: zs = dd->Zs;
1458: xm = dd->Xe-dd->Xs;
1459: ym = dd->Ye-dd->Ys;
1460: zm = dd->Ze-dd->Zs;
1461: } else {
1462: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1463: if (dd->arrayin[i]) {
1464: *iptr = dd->arrayin[i];
1465: iarray_start = (char*)dd->startin[i];
1466: dd->arrayin[i] = NULL;
1467: dd->startin[i] = NULL;
1469: goto done;
1470: }
1471: }
1472: xs = dd->xs;
1473: ys = dd->ys;
1474: zs = dd->zs;
1475: xm = dd->xe-dd->xs;
1476: ym = dd->ye-dd->ys;
1477: zm = dd->ze-dd->zs;
1478: }
1480: switch (da->dim) {
1481: case 1: {
1482: void *ptr;
1484: PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);
1486: ptr = (void*)(iarray_start - xs*sizeof(PetscScalar));
1487: *iptr = (void*)ptr;
1488: break;
1489: }
1490: case 2: {
1491: void **ptr;
1493: PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);
1495: ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
1496: for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
1497: *iptr = (void*)ptr;
1498: break;
1499: }
1500: case 3: {
1501: void ***ptr,**bptr;
1503: PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);
1505: ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1506: bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1507: for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1508: for (i=zs; i<zs+zm; i++) {
1509: for (j=ys; j<ys+ym; j++) {
1510: ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1511: }
1512: }
1513: *iptr = (void*)ptr;
1514: break;
1515: }
1516: default:
1517: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
1518: }
1520: done:
1521: /* add arrays to the checked out list */
1522: if (ghosted) {
1523: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1524: if (!dd->arrayghostedout[i]) {
1525: dd->arrayghostedout[i] = *iptr;
1526: dd->startghostedout[i] = iarray_start;
1527: break;
1528: }
1529: }
1530: } else {
1531: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1532: if (!dd->arrayout[i]) {
1533: dd->arrayout[i] = *iptr;
1534: dd->startout[i] = iarray_start;
1535: break;
1536: }
1537: }
1538: }
1539: return(0);
1540: }
1542: /*@C
1543: DMDARestoreArray - Restores an array of derivative types for a DMDA
1545: Input Parameter:
1546: + da - information about my local patch
1547: . ghosted - do you want arrays for the ghosted or nonghosted patch
1548: - vptr - array data structured to be passed to ad_FormFunctionLocal()
1550: Level: advanced
1552: .seealso: DMDAGetArray()
1554: @*/
1555: PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1556: {
1557: PetscInt i;
1558: void **iptr = (void**)vptr,*iarray_start = 0;
1559: DM_DA *dd = (DM_DA*)da->data;
1563: if (ghosted) {
1564: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1565: if (dd->arrayghostedout[i] == *iptr) {
1566: iarray_start = dd->startghostedout[i];
1567: dd->arrayghostedout[i] = NULL;
1568: dd->startghostedout[i] = NULL;
1569: break;
1570: }
1571: }
1572: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1573: if (!dd->arrayghostedin[i]) {
1574: dd->arrayghostedin[i] = *iptr;
1575: dd->startghostedin[i] = iarray_start;
1576: break;
1577: }
1578: }
1579: } else {
1580: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1581: if (dd->arrayout[i] == *iptr) {
1582: iarray_start = dd->startout[i];
1583: dd->arrayout[i] = NULL;
1584: dd->startout[i] = NULL;
1585: break;
1586: }
1587: }
1588: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1589: if (!dd->arrayin[i]) {
1590: dd->arrayin[i] = *iptr;
1591: dd->startin[i] = iarray_start;
1592: break;
1593: }
1594: }
1595: }
1596: return(0);
1597: }