Actual source code: dalocal.c
petsc-3.9.4 2018-09-11
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, numPoints, dim, dimEmbed, spDim, totDim, numValues, cStart, cEnd, f, c, v, d;
1199: PetscErrorCode ierr;
1202: DMGetDS(dm, &prob);
1203: PetscDSGetTotalDimension(prob, &totDim);
1204: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1205: PetscFEGetQuadrature(fe, &q);
1206: DMGetDefaultSection(dm, §ion);
1207: PetscSectionGetNumFields(section, &numFields);
1208: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1209: DMGetCoordinateDim(dm, &dimEmbed);
1210: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1211: DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
1212: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
1213: DMGetWorkArray(dm, numValues, MPIU_SCALAR, &values);
1214: PetscQuadratureGetData(q, NULL, NULL, &numPoints, NULL, NULL);
1215: for (c = cStart; c < cEnd; ++c) {
1216: PetscFECellGeom geom;
1218: DMDAComputeCellGeometryFEM(dm, c, q, geom.v0, geom.J, NULL, &geom.detJ);
1219: geom.dim = dim;
1220: geom.dimEmbed = dimEmbed;
1221: for (f = 0, v = 0; f < numFields; ++f) {
1222: void * const ctx = ctxs ? ctxs[f] : NULL;
1224: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1225: PetscFEGetNumComponents(fe, &Nc);
1226: PetscFEGetDualSpace(fe, &sp);
1227: PetscDualSpaceGetDimension(sp, &spDim);
1228: for (d = 0; d < spDim; ++d) {
1229: PetscDualSpaceApply(sp, d, time, &geom, Nc, funcs[f], ctx, &values[v]);
1230: }
1231: }
1232: DMDAVecSetClosure(dm, section, localX, c, values, mode);
1233: }
1234: DMRestoreWorkArray(dm, numValues, MPIU_SCALAR, &values);
1235: return(0);
1236: }
1238: PetscErrorCode DMComputeL2Diff_DA(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1239: {
1240: const PetscInt debug = 0;
1241: PetscDS prob;
1242: PetscFE fe;
1243: PetscQuadrature quad;
1244: PetscSection section;
1245: Vec localX;
1246: PetscScalar *funcVal;
1247: PetscReal *coords, *v0, *J, *invJ, detJ;
1248: PetscReal localDiff = 0.0;
1249: PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1250: PetscErrorCode ierr;
1253: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1254: DMGetDS(dm, &prob);
1255: PetscDSGetTotalComponents(prob, &Nc);
1256: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1257: PetscFEGetQuadrature(fe, &quad);
1258: DMGetDefaultSection(dm, §ion);
1259: PetscSectionGetNumFields(section, &numFields);
1260: DMGetLocalVector(dm, &localX);
1261: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1262: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1263: /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1264: PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1265: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1266: for (c = cStart; c < cEnd; ++c) {
1267: PetscScalar *x = NULL;
1268: PetscReal elemDiff = 0.0;
1270: DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);
1271: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1272: DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);
1274: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1275: void * const ctx = ctxs ? ctxs[field] : NULL;
1276: const PetscReal *quadPoints, *quadWeights;
1277: PetscReal *basis;
1278: PetscInt qNc, Nq, Nb, Nc, q, d, e, fc, f;
1280: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1281: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1282: PetscFEGetDimension(fe, &Nb);
1283: PetscFEGetNumComponents(fe, &Nc);
1284: if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Quadrature components %D != %D field components\n", qNc, Nc);
1285: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
1286: if (debug) {
1287: char title[1024];
1288: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1289: DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1290: }
1291: for (q = 0; q < Nq; ++q) {
1292: for (d = 0; d < dim; d++) {
1293: coords[d] = v0[d];
1294: for (e = 0; e < dim; e++) {
1295: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1296: }
1297: }
1298: (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);
1299: for (fc = 0; fc < Nc; ++fc) {
1300: PetscScalar interpolant = 0.0;
1302: for (f = 0; f < Nb; ++f) {
1303: interpolant += x[fieldOffset+f]*basis[(q*Nb+f)*Nc+fc];
1304: }
1305: 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);}
1306: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ;
1307: }
1308: }
1309: comp += Nc;
1310: fieldOffset += Nb;
1311: }
1312: DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1313: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
1314: localDiff += elemDiff;
1315: }
1316: PetscFree5(funcVal,coords,v0,J,invJ);
1317: DMRestoreLocalVector(dm, &localX);
1318: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1319: *diff = PetscSqrtReal(*diff);
1320: return(0);
1321: }
1323: 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)
1324: {
1325: const PetscInt debug = 0;
1326: PetscDS prob;
1327: PetscFE fe;
1328: PetscQuadrature quad;
1329: PetscSection section;
1330: Vec localX;
1331: PetscScalar *funcVal, *interpolantVec;
1332: PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
1333: PetscReal localDiff = 0.0;
1334: PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1335: PetscErrorCode ierr;
1338: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1339: DMGetDS(dm, &prob);
1340: PetscDSGetTotalComponents(prob, &Nc);
1341: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1342: PetscFEGetQuadrature(fe, &quad);
1343: DMGetDefaultSection(dm, §ion);
1344: PetscSectionGetNumFields(section, &numFields);
1345: DMGetLocalVector(dm, &localX);
1346: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1347: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1348: /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1349: PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
1350: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1351: for (c = cStart; c < cEnd; ++c) {
1352: PetscScalar *x = NULL;
1353: PetscReal elemDiff = 0.0;
1355: DMDAComputeCellGeometryFEM(dm, c, quad, v0, J, invJ, &detJ);
1356: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1357: DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);
1359: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1360: void * const ctx = ctxs ? ctxs[field] : NULL;
1361: const PetscReal *quadPoints, *quadWeights;
1362: PetscReal *basisDer;
1363: PetscInt qNc, Nq, Nb, Nc, q, d, e, fc, f, g;
1365: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1366: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1367: PetscFEGetDimension(fe, &Nb);
1368: PetscFEGetNumComponents(fe, &Nc);
1369: if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Quadrature components %D != %D field components\n", qNc, Nc);
1370: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1371: if (debug) {
1372: char title[1024];
1373: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1374: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1375: }
1376: for (q = 0; q < Nq; ++q) {
1377: for (d = 0; d < dim; d++) {
1378: coords[d] = v0[d];
1379: for (e = 0; e < dim; e++) {
1380: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1381: }
1382: }
1383: (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);
1384: for (fc = 0; fc < Nc; ++fc) {
1385: PetscScalar interpolant = 0.0;
1387: for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
1388: for (f = 0; f < Nb; ++f) {
1389: for (d = 0; d < dim; ++d) {
1390: realSpaceDer[d] = 0.0;
1391: for (g = 0; g < dim; ++g) {
1392: realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb+f)*dim+g];
1393: }
1394: interpolantVec[d] += x[fieldOffset+f]*realSpaceDer[d];
1395: }
1396: }
1397: for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
1398: 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);}
1399: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q*Nc+c]*detJ;
1400: }
1401: }
1402: comp += Nc;
1403: fieldOffset += Nb*Nc;
1404: }
1405: DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1406: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
1407: localDiff += elemDiff;
1408: }
1409: PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
1410: DMRestoreLocalVector(dm, &localX);
1411: MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1412: *diff = PetscSqrtReal(*diff);
1413: return(0);
1414: }
1416: /* ------------------------------------------------------------------- */
1418: /*@C
1419: DMDAGetArray - Gets a work array for a DMDA
1421: Input Parameter:
1422: + da - information about my local patch
1423: - ghosted - do you want arrays for the ghosted or nonghosted patch
1425: Output Parameters:
1426: . vptr - array data structured
1428: Note: The vector values are NOT initialized and may have garbage in them, so you may need
1429: to zero them.
1431: Level: advanced
1433: .seealso: DMDARestoreArray()
1435: @*/
1436: PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
1437: {
1439: PetscInt j,i,xs,ys,xm,ym,zs,zm;
1440: char *iarray_start;
1441: void **iptr = (void**)vptr;
1442: DM_DA *dd = (DM_DA*)da->data;
1446: if (ghosted) {
1447: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1448: if (dd->arrayghostedin[i]) {
1449: *iptr = dd->arrayghostedin[i];
1450: iarray_start = (char*)dd->startghostedin[i];
1451: dd->arrayghostedin[i] = NULL;
1452: dd->startghostedin[i] = NULL;
1454: goto done;
1455: }
1456: }
1457: xs = dd->Xs;
1458: ys = dd->Ys;
1459: zs = dd->Zs;
1460: xm = dd->Xe-dd->Xs;
1461: ym = dd->Ye-dd->Ys;
1462: zm = dd->Ze-dd->Zs;
1463: } else {
1464: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1465: if (dd->arrayin[i]) {
1466: *iptr = dd->arrayin[i];
1467: iarray_start = (char*)dd->startin[i];
1468: dd->arrayin[i] = NULL;
1469: dd->startin[i] = NULL;
1471: goto done;
1472: }
1473: }
1474: xs = dd->xs;
1475: ys = dd->ys;
1476: zs = dd->zs;
1477: xm = dd->xe-dd->xs;
1478: ym = dd->ye-dd->ys;
1479: zm = dd->ze-dd->zs;
1480: }
1482: switch (da->dim) {
1483: case 1: {
1484: void *ptr;
1486: PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);
1488: ptr = (void*)(iarray_start - xs*sizeof(PetscScalar));
1489: *iptr = (void*)ptr;
1490: break;
1491: }
1492: case 2: {
1493: void **ptr;
1495: PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);
1497: ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
1498: for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
1499: *iptr = (void*)ptr;
1500: break;
1501: }
1502: case 3: {
1503: void ***ptr,**bptr;
1505: PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);
1507: ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1508: bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1509: for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1510: for (i=zs; i<zs+zm; i++) {
1511: for (j=ys; j<ys+ym; j++) {
1512: ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1513: }
1514: }
1515: *iptr = (void*)ptr;
1516: break;
1517: }
1518: default:
1519: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
1520: }
1522: done:
1523: /* add arrays to the checked out list */
1524: if (ghosted) {
1525: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1526: if (!dd->arrayghostedout[i]) {
1527: dd->arrayghostedout[i] = *iptr;
1528: dd->startghostedout[i] = iarray_start;
1529: break;
1530: }
1531: }
1532: } else {
1533: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1534: if (!dd->arrayout[i]) {
1535: dd->arrayout[i] = *iptr;
1536: dd->startout[i] = iarray_start;
1537: break;
1538: }
1539: }
1540: }
1541: return(0);
1542: }
1544: /*@C
1545: DMDARestoreArray - Restores an array of derivative types for a DMDA
1547: Input Parameter:
1548: + da - information about my local patch
1549: . ghosted - do you want arrays for the ghosted or nonghosted patch
1550: - vptr - array data structured to be passed to ad_FormFunctionLocal()
1552: Level: advanced
1554: .seealso: DMDAGetArray()
1556: @*/
1557: PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1558: {
1559: PetscInt i;
1560: void **iptr = (void**)vptr,*iarray_start = 0;
1561: DM_DA *dd = (DM_DA*)da->data;
1565: if (ghosted) {
1566: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1567: if (dd->arrayghostedout[i] == *iptr) {
1568: iarray_start = dd->startghostedout[i];
1569: dd->arrayghostedout[i] = NULL;
1570: dd->startghostedout[i] = NULL;
1571: break;
1572: }
1573: }
1574: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1575: if (!dd->arrayghostedin[i]) {
1576: dd->arrayghostedin[i] = *iptr;
1577: dd->startghostedin[i] = iarray_start;
1578: break;
1579: }
1580: }
1581: } else {
1582: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1583: if (dd->arrayout[i] == *iptr) {
1584: iarray_start = dd->startout[i];
1585: dd->arrayout[i] = NULL;
1586: dd->startout[i] = NULL;
1587: break;
1588: }
1589: }
1590: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1591: if (!dd->arrayin[i]) {
1592: dd->arrayin[i] = *iptr;
1593: dd->startin[i] = iarray_start;
1594: break;
1595: }
1596: }
1597: }
1598: return(0);
1599: }