Actual source code: dalocal.c
petsc-3.5.4 2015-05-23
2: /*
3: Code for manipulating distributed regular arrays in parallel.
4: */
6: #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/
7: #include <petscbt.h>
8: #include <petscsf.h>
9: #include <petscds.h>
10: #include <petscfe.h>
12: /*
13: This allows the DMDA vectors to properly tell MATLAB their dimensions
14: */
15: #if defined(PETSC_HAVE_MATLAB_ENGINE)
16: #include <engine.h> /* MATLAB include file */
17: #include <mex.h> /* MATLAB include file */
20: static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
21: {
23: PetscInt n,m;
24: Vec vec = (Vec)obj;
25: PetscScalar *array;
26: mxArray *mat;
27: DM da;
30: VecGetDM(vec, &da);
31: if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA");
32: DMDAGetGhostCorners(da,0,0,0,&m,&n,0);
34: VecGetArray(vec,&array);
35: #if !defined(PETSC_USE_COMPLEX)
36: mat = mxCreateDoubleMatrix(m,n,mxREAL);
37: #else
38: mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
39: #endif
40: PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));
41: PetscObjectName(obj);
42: engPutVariable((Engine*)mengine,obj->name,mat);
44: VecRestoreArray(vec,&array);
45: return(0);
46: }
47: #endif
52: PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g)
53: {
55: DM_DA *dd = (DM_DA*)da->data;
60: if (da->defaultSection) {
61: DMCreateLocalVector_Section_Private(da,g);
62: } else {
63: VecCreate(PETSC_COMM_SELF,g);
64: VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);
65: VecSetBlockSize(*g,dd->w);
66: VecSetType(*g,da->vectype);
67: VecSetDM(*g, da);
68: #if defined(PETSC_HAVE_MATLAB_ENGINE)
69: if (dd->w == 1 && dd->dim == 2) {
70: PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);
71: }
72: #endif
73: }
74: return(0);
75: }
79: /*@
80: DMDAGetNumCells - Get the number of cells in the local piece of the DMDA. This includes ghost cells.
82: Input Parameter:
83: . dm - The DM object
85: Output Parameters:
86: + numCellsX - The number of local cells in the x-direction
87: . numCellsY - The number of local cells in the y-direction
88: . numCellsZ - The number of local cells in the z-direction
89: - numCells - The number of local cells
91: Level: developer
93: .seealso: DMDAGetCellPoint()
94: @*/
95: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
96: {
97: DM_DA *da = (DM_DA*) dm->data;
98: const PetscInt dim = da->dim;
99: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
100: const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
104: if (numCellsX) {
106: *numCellsX = mx;
107: }
108: if (numCellsY) {
110: *numCellsY = my;
111: }
112: if (numCellsZ) {
114: *numCellsZ = mz;
115: }
116: if (numCells) {
118: *numCells = nC;
119: }
120: return(0);
121: }
125: /*@
126: DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the DMDA
128: Input Parameters:
129: + dm - The DM object
130: - i,j,k - The global indices for the cell
132: Output Parameters:
133: . point - The local DM point
135: Level: developer
137: .seealso: DMDAGetNumCells()
138: @*/
139: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
140: {
141: DM_DA *da = (DM_DA*) dm->data;
142: const PetscInt dim = da->dim;
143: DMDALocalInfo info;
149: DMDAGetLocalInfo(dm, &info);
150: if (dim > 0) {if ((i < info.gxs) || (i >= info.gxs+info.gxm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "X index %d not in [%d, %d)", i, info.gxs, info.gxs+info.gxm);}
151: if (dim > 1) {if ((i < info.gys) || (i >= info.gys+info.gym)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %d not in [%d, %d)", i, info.gys, info.gys+info.gym);}
152: if (dim > 2) {if ((i < info.gzs) || (i >= info.gzs+info.gzm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %d not in [%d, %d)", i, info.gzs, info.gzs+info.gzm);}
153: *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0);
154: return(0);
155: }
159: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
160: {
161: DM_DA *da = (DM_DA*) dm->data;
162: const PetscInt dim = da->dim;
163: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
164: const PetscInt nVx = mx+1;
165: const PetscInt nVy = dim > 1 ? (my+1) : 1;
166: const PetscInt nVz = dim > 2 ? (mz+1) : 1;
167: const PetscInt nV = nVx*nVy*nVz;
170: if (numVerticesX) {
172: *numVerticesX = nVx;
173: }
174: if (numVerticesY) {
176: *numVerticesY = nVy;
177: }
178: if (numVerticesZ) {
180: *numVerticesZ = nVz;
181: }
182: if (numVertices) {
184: *numVertices = nV;
185: }
186: return(0);
187: }
191: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
192: {
193: DM_DA *da = (DM_DA*) dm->data;
194: const PetscInt dim = da->dim;
195: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
196: const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
197: const PetscInt nXF = (mx+1)*nxF;
198: const PetscInt nyF = mx*(dim > 2 ? mz : 1);
199: const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0;
200: const PetscInt nzF = mx*(dim > 1 ? my : 0);
201: const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0;
204: if (numXFacesX) {
206: *numXFacesX = nxF;
207: }
208: if (numXFaces) {
210: *numXFaces = nXF;
211: }
212: if (numYFacesY) {
214: *numYFacesY = nyF;
215: }
216: if (numYFaces) {
218: *numYFaces = nYF;
219: }
220: if (numZFacesZ) {
222: *numZFacesZ = nzF;
223: }
224: if (numZFaces) {
226: *numZFaces = nZF;
227: }
228: return(0);
229: }
233: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
234: {
235: DM_DA *da = (DM_DA*) dm->data;
236: const PetscInt dim = da->dim;
237: PetscInt nC, nV, nXF, nYF, nZF;
243: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
244: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
245: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
246: if (height == 0) {
247: /* Cells */
248: if (pStart) *pStart = 0;
249: if (pEnd) *pEnd = nC;
250: } else if (height == 1) {
251: /* Faces */
252: if (pStart) *pStart = nC+nV;
253: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
254: } else if (height == dim) {
255: /* Vertices */
256: if (pStart) *pStart = nC;
257: if (pEnd) *pEnd = nC+nV;
258: } else if (height < 0) {
259: /* All points */
260: if (pStart) *pStart = 0;
261: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
262: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
263: return(0);
264: }
268: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
269: {
270: DM_DA *da = (DM_DA*) dm->data;
271: const PetscInt dim = da->dim;
272: PetscInt nC, nV, nXF, nYF, nZF;
278: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
279: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
280: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
281: if (depth == dim) {
282: /* Cells */
283: if (pStart) *pStart = 0;
284: if (pEnd) *pEnd = nC;
285: } else if (depth == dim-1) {
286: /* Faces */
287: if (pStart) *pStart = nC+nV;
288: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
289: } else if (depth == 0) {
290: /* Vertices */
291: if (pStart) *pStart = nC;
292: if (pEnd) *pEnd = nC+nV;
293: } else if (depth < 0) {
294: /* All points */
295: if (pStart) *pStart = 0;
296: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
297: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
298: return(0);
299: }
303: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
304: {
305: DM_DA *da = (DM_DA*) dm->data;
306: const PetscInt dim = da->dim;
307: PetscInt nC, nV, nXF, nYF, nZF;
311: *coneSize = 0;
312: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
313: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
314: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
315: switch (dim) {
316: case 2:
317: if (p >= 0) {
318: if (p < nC) {
319: *coneSize = 4;
320: } else if (p < nC+nV) {
321: *coneSize = 0;
322: } else if (p < nC+nV+nXF+nYF+nZF) {
323: *coneSize = 2;
324: } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
325: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
326: break;
327: case 3:
328: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
329: break;
330: }
331: return(0);
332: }
336: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
337: {
338: DM_DA *da = (DM_DA*) dm->data;
339: const PetscInt dim = da->dim;
340: PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
344: if (!cone) {DMGetWorkArray(dm, 6, PETSC_INT, cone);}
345: DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);
346: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
347: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
348: switch (dim) {
349: case 2:
350: if (p >= 0) {
351: if (p < nC) {
352: const PetscInt cy = p / nCx;
353: const PetscInt cx = p % nCx;
355: (*cone)[0] = cy*nxF + cx + nC+nV;
356: (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
357: (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
358: (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
359: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
360: } else if (p < nC+nV) {
361: } else if (p < nC+nV+nXF) {
362: const PetscInt fy = (p - nC+nV) / nxF;
363: const PetscInt fx = (p - nC+nV) % nxF;
365: (*cone)[0] = fy*nVx + fx + nC;
366: (*cone)[1] = fy*nVx + fx + 1 + nC;
367: } else if (p < nC+nV+nXF+nYF) {
368: const PetscInt fx = (p - nC+nV+nXF) / nyF;
369: const PetscInt fy = (p - nC+nV+nXF) % nyF;
371: (*cone)[0] = fy*nVx + fx + nC;
372: (*cone)[1] = fy*nVx + fx + nVx + nC;
373: } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
374: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
375: break;
376: case 3:
377: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
378: break;
379: }
380: return(0);
381: }
385: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
386: {
390: DMGetWorkArray(dm, 6, PETSC_INT, cone);
391: return(0);
392: }
396: /*@C
397: DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with
398: different numbers of dofs on vertices, cells, and faces in each direction.
400: Input Parameters:
401: + dm- The DMDA
402: . numFields - The number of fields
403: . numComp - The number of components in each field
404: . numDof - The number of dofs per dimension for each field
405: . numFaceDof - The number of dofs per face for each field and direction, or NULL
407: Level: developer
409: Note:
410: The default DMDA numbering is as follows:
412: - Cells: [0, nC)
413: - Vertices: [nC, nC+nV)
414: - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir
415: - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir
416: - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir
418: We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment.
419: @*/
420: PetscErrorCode DMDACreateSection(DM dm, const PetscInt numComp[], const PetscInt numDof[], const PetscInt numFaceDof[], PetscSection *s)
421: {
422: DM_DA *da = (DM_DA*) dm->data;
423: PetscSection section;
424: const PetscInt dim = da->dim;
425: PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0};
426: PetscBT isLeaf;
427: PetscSF sf;
428: PetscMPIInt rank;
429: const PetscMPIInt *neighbors;
430: PetscInt *localPoints;
431: PetscSFNode *remotePoints;
432: PetscInt nleaves = 0, nleavesCheck = 0, nL = 0;
433: PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF;
434: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd;
435: PetscInt f, v, c, xf, yf, zf, xn, yn, zn;
436: PetscErrorCode ierr;
441: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
442: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
443: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
444: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
445: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
446: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
447: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
448: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
449: xfStart = vEnd; xfEnd = xfStart+nXF;
450: yfStart = xfEnd; yfEnd = yfStart+nYF;
451: zfStart = yfEnd; zfEnd = zfStart+nZF;
452: if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd);
453: /* Create local section */
454: DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);
455: for (f = 0; f < numFields; ++f) {
456: numVertexTotDof += numDof[f*(dim+1)+0];
457: numCellTotDof += numDof[f*(dim+1)+dim];
458: numFaceTotDof[0] += dim > 0 ? (numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]) : 0;
459: numFaceTotDof[1] += dim > 1 ? (numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]) : 0;
460: numFaceTotDof[2] += dim > 2 ? (numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]) : 0;
461: }
462: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
463: if (numFields > 0) {
464: PetscSectionSetNumFields(section, numFields);
465: for (f = 0; f < numFields; ++f) {
466: const char *name;
468: DMDAGetFieldName(dm, f, &name);
469: PetscSectionSetFieldName(section, f, name ? name : "Field");
470: if (numComp) {
471: PetscSectionSetFieldComponents(section, f, numComp[f]);
472: }
473: }
474: }
475: PetscSectionSetChart(section, pStart, pEnd);
476: for (v = vStart; v < vEnd; ++v) {
477: for (f = 0; f < numFields; ++f) {
478: PetscSectionSetFieldDof(section, v, f, numDof[f*(dim+1)+0]);
479: }
480: PetscSectionSetDof(section, v, numVertexTotDof);
481: }
482: for (xf = xfStart; xf < xfEnd; ++xf) {
483: for (f = 0; f < numFields; ++f) {
484: PetscSectionSetFieldDof(section, xf, f, numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]);
485: }
486: PetscSectionSetDof(section, xf, numFaceTotDof[0]);
487: }
488: for (yf = yfStart; yf < yfEnd; ++yf) {
489: for (f = 0; f < numFields; ++f) {
490: PetscSectionSetFieldDof(section, yf, f, numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]);
491: }
492: PetscSectionSetDof(section, yf, numFaceTotDof[1]);
493: }
494: for (zf = zfStart; zf < zfEnd; ++zf) {
495: for (f = 0; f < numFields; ++f) {
496: PetscSectionSetFieldDof(section, zf, f, numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]);
497: }
498: PetscSectionSetDof(section, zf, numFaceTotDof[2]);
499: }
500: for (c = cStart; c < cEnd; ++c) {
501: for (f = 0; f < numFields; ++f) {
502: PetscSectionSetFieldDof(section, c, f, numDof[f*(dim+1)+dim]);
503: }
504: PetscSectionSetDof(section, c, numCellTotDof);
505: }
506: PetscSectionSetUp(section);
507: /* Create mesh point SF */
508: PetscBTCreate(pEnd-pStart, &isLeaf);
509: DMDAGetNeighbors(dm, &neighbors);
510: for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
511: for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
512: for (xn = 0; xn < 3; ++xn) {
513: const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
514: const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
515: PetscInt xv, yv, zv;
517: if (neighbor >= 0 && neighbor < rank) {
518: if (xp < 0) { /* left */
519: if (yp < 0) { /* bottom */
520: if (zp < 0) { /* back */
521: const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* left bottom back vertex */
522: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
523: } else if (zp > 0) { /* front */
524: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* left bottom front vertex */
525: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
526: } else {
527: for (zv = 0; zv < nVz; ++zv) {
528: const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; /* left bottom vertices */
529: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
530: }
531: }
532: } else if (yp > 0) { /* top */
533: if (zp < 0) { /* back */
534: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* left top back vertex */
535: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
536: } else if (zp > 0) { /* front */
537: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* left top front vertex */
538: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
539: } else {
540: for (zv = 0; zv < nVz; ++zv) {
541: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* left top vertices */
542: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
543: }
544: }
545: } else {
546: if (zp < 0) { /* back */
547: for (yv = 0; yv < nVy; ++yv) {
548: const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* left back vertices */
549: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
550: }
551: } else if (zp > 0) { /* front */
552: for (yv = 0; yv < nVy; ++yv) {
553: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* left front vertices */
554: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
555: }
556: } else {
557: for (zv = 0; zv < nVz; ++zv) {
558: for (yv = 0; yv < nVy; ++yv) {
559: const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; /* left vertices */
560: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
561: }
562: }
563: #if 0
564: for (xf = 0; xf < nxF; ++xf) {
565: /* THIS IS WRONG */
566: const PetscInt localFace = 0 + nC+nV; /* left faces */
567: if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
568: }
569: #endif
570: }
571: }
572: } else if (xp > 0) { /* right */
573: if (yp < 0) { /* bottom */
574: if (zp < 0) { /* back */
575: const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* right bottom back vertex */
576: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
577: } else if (zp > 0) { /* front */
578: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* right bottom front vertex */
579: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
580: } else {
581: for (zv = 0; zv < nVz; ++zv) {
582: const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* right bottom vertices */
583: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
584: }
585: }
586: } else if (yp > 0) { /* top */
587: if (zp < 0) { /* back */
588: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
589: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
590: } else if (zp > 0) { /* front */
591: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
592: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
593: } else {
594: for (zv = 0; zv < nVz; ++zv) {
595: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
596: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
597: }
598: }
599: } else {
600: if (zp < 0) { /* back */
601: for (yv = 0; yv < nVy; ++yv) {
602: const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
603: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
604: }
605: } else if (zp > 0) { /* front */
606: for (yv = 0; yv < nVy; ++yv) {
607: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
608: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
609: }
610: } else {
611: for (zv = 0; zv < nVz; ++zv) {
612: for (yv = 0; yv < nVy; ++yv) {
613: const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
614: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
615: }
616: }
617: #if 0
618: for (xf = 0; xf < nxF; ++xf) {
619: /* THIS IS WRONG */
620: const PetscInt localFace = 0 + nC+nV; /* right faces */
621: if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
622: }
623: #endif
624: }
625: }
626: } else {
627: if (yp < 0) { /* bottom */
628: if (zp < 0) { /* back */
629: for (xv = 0; xv < nVx; ++xv) {
630: const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; /* bottom back vertices */
631: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
632: }
633: } else if (zp > 0) { /* front */
634: for (xv = 0; xv < nVx; ++xv) {
635: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* bottom front vertices */
636: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
637: }
638: } else {
639: for (zv = 0; zv < nVz; ++zv) {
640: for (xv = 0; xv < nVx; ++xv) {
641: const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; /* bottom vertices */
642: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
643: }
644: }
645: #if 0
646: for (yf = 0; yf < nyF; ++yf) {
647: /* THIS IS WRONG */
648: const PetscInt localFace = 0 + nC+nV; /* bottom faces */
649: if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves;
650: }
651: #endif
652: }
653: } else if (yp > 0) { /* top */
654: if (zp < 0) { /* back */
655: for (xv = 0; xv < nVx; ++xv) {
656: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
657: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
658: }
659: } else if (zp > 0) { /* front */
660: for (xv = 0; xv < nVx; ++xv) {
661: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
662: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
663: }
664: } else {
665: for (zv = 0; zv < nVz; ++zv) {
666: for (xv = 0; xv < nVx; ++xv) {
667: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
668: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
669: }
670: }
671: #if 0
672: for (yf = 0; yf < nyF; ++yf) {
673: /* THIS IS WRONG */
674: const PetscInt localFace = 0 + nC+nV; /* top faces */
675: if (!PetscBTLookupSet(isLeaf, localVFace)) ++nleaves;
676: }
677: #endif
678: }
679: } else {
680: if (zp < 0) { /* back */
681: for (yv = 0; yv < nVy; ++yv) {
682: for (xv = 0; xv < nVx; ++xv) {
683: const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; /* back vertices */
684: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
685: }
686: }
687: #if 0
688: for (zf = 0; zf < nzF; ++zf) {
689: /* THIS IS WRONG */
690: const PetscInt localFace = 0 + nC+nV; /* back faces */
691: if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
692: }
693: #endif
694: } else if (zp > 0) { /* front */
695: for (yv = 0; yv < nVy; ++yv) {
696: for (xv = 0; xv < nVx; ++xv) {
697: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
698: if (!PetscBTLookupSet(isLeaf, localVertex)) ++nleaves;
699: }
700: }
701: #if 0
702: for (zf = 0; zf < nzF; ++zf) {
703: /* THIS IS WRONG */
704: const PetscInt localFace = 0 + nC+nV; /* front faces */
705: if (!PetscBTLookupSet(isLeaf, localFace)) ++nleaves;
706: }
707: #endif
708: } else {
709: /* Nothing is shared from the interior */
710: }
711: }
712: }
713: }
714: }
715: }
716: }
717: PetscBTMemzero(pEnd-pStart, isLeaf);
718: PetscMalloc2(nleaves,&localPoints,nleaves,&remotePoints);
719: for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
720: for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
721: for (xn = 0; xn < 3; ++xn) {
722: const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
723: const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
724: PetscInt xv, yv, zv;
726: if (neighbor >= 0 && neighbor < rank) {
727: if (xp < 0) { /* left */
728: if (yp < 0) { /* bottom */
729: if (zp < 0) { /* back */
730: const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* left bottom back vertex */
731: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
733: if (!PetscBTLookupSet(isLeaf, localVertex)) {
734: localPoints[nL] = localVertex;
735: remotePoints[nL].rank = neighbor;
736: remotePoints[nL].index = remoteVertex;
737: ++nL;
738: }
739: } else if (zp > 0) { /* front */
740: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* left bottom front vertex */
741: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
743: if (!PetscBTLookupSet(isLeaf, localVertex)) {
744: localPoints[nL] = localVertex;
745: remotePoints[nL].rank = neighbor;
746: remotePoints[nL].index = remoteVertex;
747: ++nL;
748: }
749: } else {
750: for (zv = 0; zv < nVz; ++zv) {
751: const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; /* left bottom vertices */
752: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
754: if (!PetscBTLookupSet(isLeaf, localVertex)) {
755: localPoints[nL] = localVertex;
756: remotePoints[nL].rank = neighbor;
757: remotePoints[nL].index = remoteVertex;
758: ++nL;
759: }
760: }
761: }
762: } else if (yp > 0) { /* top */
763: if (zp < 0) { /* back */
764: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* left top back vertex */
765: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
767: if (!PetscBTLookupSet(isLeaf, localVertex)) {
768: localPoints[nL] = localVertex;
769: remotePoints[nL].rank = neighbor;
770: remotePoints[nL].index = remoteVertex;
771: ++nL;
772: }
773: } else if (zp > 0) { /* front */
774: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* left top front vertex */
775: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
777: if (!PetscBTLookupSet(isLeaf, localVertex)) {
778: localPoints[nL] = localVertex;
779: remotePoints[nL].rank = neighbor;
780: remotePoints[nL].index = remoteVertex;
781: ++nL;
782: }
783: } else {
784: for (zv = 0; zv < nVz; ++zv) {
785: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* left top vertices */
786: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
788: if (!PetscBTLookupSet(isLeaf, localVertex)) {
789: localPoints[nL] = localVertex;
790: remotePoints[nL].rank = neighbor;
791: remotePoints[nL].index = remoteVertex;
792: ++nL;
793: }
794: }
795: }
796: } else {
797: if (zp < 0) { /* back */
798: for (yv = 0; yv < nVy; ++yv) {
799: const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* left back vertices */
800: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
802: if (!PetscBTLookupSet(isLeaf, localVertex)) {
803: localPoints[nL] = localVertex;
804: remotePoints[nL].rank = neighbor;
805: remotePoints[nL].index = remoteVertex;
806: ++nL;
807: }
808: }
809: } else if (zp > 0) { /* front */
810: for (yv = 0; yv < nVy; ++yv) {
811: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* left front vertices */
812: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
814: if (!PetscBTLookupSet(isLeaf, localVertex)) {
815: localPoints[nL] = localVertex;
816: remotePoints[nL].rank = neighbor;
817: remotePoints[nL].index = remoteVertex;
818: ++nL;
819: }
820: }
821: } else {
822: for (zv = 0; zv < nVz; ++zv) {
823: for (yv = 0; yv < nVy; ++yv) {
824: const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; /* left vertices */
825: const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
827: if (!PetscBTLookupSet(isLeaf, localVertex)) {
828: localPoints[nL] = localVertex;
829: remotePoints[nL].rank = neighbor;
830: remotePoints[nL].index = remoteVertex;
831: ++nL;
832: }
833: }
834: }
835: #if 0
836: for (xf = 0; xf < nxF; ++xf) {
837: /* THIS IS WRONG */
838: const PetscInt localFace = 0 + nC+nV; /* left faces */
839: const PetscInt remoteFace = 0 + nC+nV;
841: if (!PetscBTLookupSet(isLeaf, localFace)) {
842: localPoints[nL] = localFace;
843: remotePoints[nL].rank = neighbor;
844: remotePoints[nL].index = remoteFace;
845: }
846: }
847: #endif
848: }
849: }
850: } else if (xp > 0) { /* right */
851: if (yp < 0) { /* bottom */
852: if (zp < 0) { /* back */
853: const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* right bottom back vertex */
854: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
856: if (!PetscBTLookupSet(isLeaf, localVertex)) {
857: localPoints[nL] = localVertex;
858: remotePoints[nL].rank = neighbor;
859: remotePoints[nL].index = remoteVertex;
860: ++nL;
861: }
862: } else if (zp > 0) { /* front */
863: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* right bottom front vertex */
864: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
866: if (!PetscBTLookupSet(isLeaf, localVertex)) {
867: localPoints[nL] = localVertex;
868: remotePoints[nL].rank = neighbor;
869: remotePoints[nL].index = remoteVertex;
870: ++nL;
871: }
872: } else {
873: nleavesCheck += nVz;
874: for (zv = 0; zv < nVz; ++zv) {
875: const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* right bottom vertices */
876: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
878: if (!PetscBTLookupSet(isLeaf, localVertex)) {
879: localPoints[nL] = localVertex;
880: remotePoints[nL].rank = neighbor;
881: remotePoints[nL].index = remoteVertex;
882: ++nL;
883: }
884: }
885: }
886: } else if (yp > 0) { /* top */
887: if (zp < 0) { /* back */
888: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top back vertex */
889: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
891: if (!PetscBTLookupSet(isLeaf, localVertex)) {
892: localPoints[nL] = localVertex;
893: remotePoints[nL].rank = neighbor;
894: remotePoints[nL].index = remoteVertex;
895: ++nL;
896: }
897: } else if (zp > 0) { /* front */
898: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top front vertex */
899: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
901: if (!PetscBTLookupSet(isLeaf, localVertex)) {
902: localPoints[nL] = localVertex;
903: remotePoints[nL].rank = neighbor;
904: remotePoints[nL].index = remoteVertex;
905: ++nL;
906: }
907: } else {
908: for (zv = 0; zv < nVz; ++zv) {
909: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* right top vertices */
910: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
912: if (!PetscBTLookupSet(isLeaf, localVertex)) {
913: localPoints[nL] = localVertex;
914: remotePoints[nL].rank = neighbor;
915: remotePoints[nL].index = remoteVertex;
916: ++nL;
917: }
918: }
919: }
920: } else {
921: if (zp < 0) { /* back */
922: for (yv = 0; yv < nVy; ++yv) {
923: const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* right back vertices */
924: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
926: if (!PetscBTLookupSet(isLeaf, localVertex)) {
927: localPoints[nL] = localVertex;
928: remotePoints[nL].rank = neighbor;
929: remotePoints[nL].index = remoteVertex;
930: ++nL;
931: }
932: }
933: } else if (zp > 0) { /* front */
934: for (yv = 0; yv < nVy; ++yv) {
935: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* right front vertices */
936: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
938: if (!PetscBTLookupSet(isLeaf, localVertex)) {
939: localPoints[nL] = localVertex;
940: remotePoints[nL].rank = neighbor;
941: remotePoints[nL].index = remoteVertex;
942: ++nL;
943: }
944: }
945: } else {
946: for (zv = 0; zv < nVz; ++zv) {
947: for (yv = 0; yv < nVy; ++yv) {
948: const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* right vertices */
949: const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
951: if (!PetscBTLookupSet(isLeaf, localVertex)) {
952: localPoints[nL] = localVertex;
953: remotePoints[nL].rank = neighbor;
954: remotePoints[nL].index = remoteVertex;
955: ++nL;
956: }
957: }
958: }
959: #if 0
960: for (xf = 0; xf < nxF; ++xf) {
961: /* THIS IS WRONG */
962: const PetscInt localFace = 0 + nC+nV; /* right faces */
963: const PetscInt remoteFace = 0 + nC+nV;
965: if (!PetscBTLookupSet(isLeaf, localFace)) {
966: localPoints[nL] = localFace;
967: remotePoints[nL].rank = neighbor;
968: remotePoints[nL].index = remoteFace;
969: ++nL;
970: }
971: }
972: #endif
973: }
974: }
975: } else {
976: if (yp < 0) { /* bottom */
977: if (zp < 0) { /* back */
978: for (xv = 0; xv < nVx; ++xv) {
979: const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; /* bottom back vertices */
980: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
982: if (!PetscBTLookupSet(isLeaf, localVertex)) {
983: localPoints[nL] = localVertex;
984: remotePoints[nL].rank = neighbor;
985: remotePoints[nL].index = remoteVertex;
986: ++nL;
987: }
988: }
989: } else if (zp > 0) { /* front */
990: for (xv = 0; xv < nVx; ++xv) {
991: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* bottom front vertices */
992: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
994: if (!PetscBTLookupSet(isLeaf, localVertex)) {
995: localPoints[nL] = localVertex;
996: remotePoints[nL].rank = neighbor;
997: remotePoints[nL].index = remoteVertex;
998: ++nL;
999: }
1000: }
1001: } else {
1002: for (zv = 0; zv < nVz; ++zv) {
1003: for (xv = 0; xv < nVx; ++xv) {
1004: const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; /* bottom vertices */
1005: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1007: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1008: localPoints[nL] = localVertex;
1009: remotePoints[nL].rank = neighbor;
1010: remotePoints[nL].index = remoteVertex;
1011: ++nL;
1012: }
1013: }
1014: }
1015: #if 0
1016: for (yf = 0; yf < nyF; ++yf) {
1017: /* THIS IS WRONG */
1018: const PetscInt localFace = 0 + nC+nV; /* bottom faces */
1019: const PetscInt remoteFace = 0 + nC+nV;
1021: if (!PetscBTLookupSet(isLeaf, localFace)) {
1022: localPoints[nL] = localFace;
1023: remotePoints[nL].rank = neighbor;
1024: remotePoints[nL].index = remoteFace;
1025: ++nL;
1026: }
1027: }
1028: #endif
1029: }
1030: } else if (yp > 0) { /* top */
1031: if (zp < 0) { /* back */
1032: for (xv = 0; xv < nVx; ++xv) {
1033: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* top back vertices */
1034: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1036: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1037: localPoints[nL] = localVertex;
1038: remotePoints[nL].rank = neighbor;
1039: remotePoints[nL].index = remoteVertex;
1040: ++nL;
1041: }
1042: }
1043: } else if (zp > 0) { /* front */
1044: for (xv = 0; xv < nVx; ++xv) {
1045: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* top front vertices */
1046: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1048: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1049: localPoints[nL] = localVertex;
1050: remotePoints[nL].rank = neighbor;
1051: remotePoints[nL].index = remoteVertex;
1052: ++nL;
1053: }
1054: }
1055: } else {
1056: for (zv = 0; zv < nVz; ++zv) {
1057: for (xv = 0; xv < nVx; ++xv) {
1058: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* top vertices */
1059: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1061: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1062: localPoints[nL] = localVertex;
1063: remotePoints[nL].rank = neighbor;
1064: remotePoints[nL].index = remoteVertex;
1065: ++nL;
1066: }
1067: }
1068: }
1069: #if 0
1070: for (yf = 0; yf < nyF; ++yf) {
1071: /* THIS IS WRONG */
1072: const PetscInt localFace = 0 + nC+nV; /* top faces */
1073: const PetscInt remoteFace = 0 + nC+nV;
1075: if (!PetscBTLookupSet(isLeaf, localFace)) {
1076: localPoints[nL] = localFace;
1077: remotePoints[nL].rank = neighbor;
1078: remotePoints[nL].index = remoteFace;
1079: ++nL;
1080: }
1081: }
1082: #endif
1083: }
1084: } else {
1085: if (zp < 0) { /* back */
1086: for (yv = 0; yv < nVy; ++yv) {
1087: for (xv = 0; xv < nVx; ++xv) {
1088: const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; /* back vertices */
1089: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1091: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1092: localPoints[nL] = localVertex;
1093: remotePoints[nL].rank = neighbor;
1094: remotePoints[nL].index = remoteVertex;
1095: ++nL;
1096: }
1097: }
1098: }
1099: #if 0
1100: for (zf = 0; zf < nzF; ++zf) {
1101: /* THIS IS WRONG */
1102: const PetscInt localFace = 0 + nC+nV; /* back faces */
1103: const PetscInt remoteFace = 0 + nC+nV;
1105: if (!PetscBTLookupSet(isLeaf, localFace)) {
1106: localPoints[nL] = localFace;
1107: remotePoints[nL].rank = neighbor;
1108: remotePoints[nL].index = remoteFace;
1109: ++nL;
1110: }
1111: }
1112: #endif
1113: } else if (zp > 0) { /* front */
1114: for (yv = 0; yv < nVy; ++yv) {
1115: for (xv = 0; xv < nVx; ++xv) {
1116: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* front vertices */
1117: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
1119: if (!PetscBTLookupSet(isLeaf, localVertex)) {
1120: localPoints[nL] = localVertex;
1121: remotePoints[nL].rank = neighbor;
1122: remotePoints[nL].index = remoteVertex;
1123: ++nL;
1124: }
1125: }
1126: }
1127: #if 0
1128: for (zf = 0; zf < nzF; ++zf) {
1129: /* THIS IS WRONG */
1130: const PetscInt localFace = 0 + nC+nV; /* front faces */
1131: const PetscInt remoteFace = 0 + nC+nV;
1133: if (!PetscBTLookupSet(isLeaf, localFace)) {
1134: localPoints[nL] = localFace;
1135: remotePoints[nL].rank = neighbor;
1136: remotePoints[nL].index = remoteFace;
1137: ++nL;
1138: }
1139: }
1140: #endif
1141: } else {
1142: /* Nothing is shared from the interior */
1143: }
1144: }
1145: }
1146: }
1147: }
1148: }
1149: }
1150: PetscBTDestroy(&isLeaf);
1151: /* Remove duplication in leaf determination */
1152: if (nleaves != nL) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck);
1153: PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);
1154: PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1155: DMSetPointSF(dm, sf);
1156: PetscSFDestroy(&sf);
1157: *s = section;
1158: return(0);
1159: }
1163: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
1164: {
1165: DM_DA *da = (DM_DA *) dm->data;
1166: Vec coordinates;
1167: PetscSection section;
1168: PetscScalar *coords;
1169: PetscReal h[3];
1170: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
1175: DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);
1176: h[0] = (xu - xl)/M;
1177: h[1] = (yu - yl)/N;
1178: h[2] = (zu - zl)/P;
1179: DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
1180: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
1181: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);
1182: PetscSectionSetNumFields(section, 1);
1183: PetscSectionSetFieldComponents(section, 0, dim);
1184: PetscSectionSetChart(section, vStart, vEnd);
1185: for (v = vStart; v < vEnd; ++v) {
1186: PetscSectionSetDof(section, v, dim);
1187: }
1188: PetscSectionSetUp(section);
1189: PetscSectionGetStorageSize(section, &size);
1190: VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
1191: VecGetArray(coordinates, &coords);
1192: for (k = 0; k < nVz; ++k) {
1193: PetscInt ind[3], d, off;
1195: ind[0] = 0;
1196: ind[1] = 0;
1197: ind[2] = k + da->zs;
1198: for (j = 0; j < nVy; ++j) {
1199: ind[1] = j + da->ys;
1200: for (i = 0; i < nVx; ++i) {
1201: const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;
1203: PetscSectionGetOffset(section, vertex, &off);
1204: ind[0] = i + da->xs;
1205: for (d = 0; d < dim; ++d) {
1206: coords[off+d] = h[d]*ind[d];
1207: }
1208: }
1209: }
1210: }
1211: VecRestoreArray(coordinates, &coords);
1212: DMSetCoordinateSection(dm, section);
1213: DMSetCoordinatesLocal(dm, coordinates);
1214: PetscSectionDestroy(§ion);
1215: VecDestroy(&coordinates);
1216: return(0);
1217: }
1221: PetscErrorCode DMDAProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
1222: {
1223: PetscDS prob;
1224: PetscFE fe;
1225: PetscDualSpace sp;
1226: PetscQuadrature q;
1227: PetscSection section;
1228: PetscScalar *values;
1229: PetscReal *v0, *J, *detJ;
1230: PetscInt numFields, numComp, numPoints, dim, spDim, totDim, numValues, cStart, cEnd, f, c, v, d;
1231: PetscErrorCode ierr;
1234: DMGetDS(dm, &prob);
1235: PetscDSGetTotalDimension(prob, &totDim);
1236: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1237: PetscFEGetQuadrature(fe, &q);
1238: DMGetDefaultSection(dm, §ion);
1239: PetscSectionGetNumFields(section, &numFields);
1240: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1241: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1242: DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
1243: if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
1244: DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
1245: PetscQuadratureGetData(q, NULL, &numPoints, NULL, NULL);
1246: PetscMalloc3(dim*numPoints,&v0,dim*dim*numPoints,&J,numPoints,&detJ);
1247: for (c = cStart; c < cEnd; ++c) {
1248: PetscCellGeometry geom;
1250: DMDAComputeCellGeometry(dm, c, q, v0, J, NULL, detJ);
1251: geom.v0 = v0;
1252: geom.J = J;
1253: geom.detJ = detJ;
1254: for (f = 0, v = 0; f < numFields; ++f) {
1255: void * const ctx = ctxs ? ctxs[f] : NULL;
1257: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1258: PetscFEGetNumComponents(fe, &numComp);
1259: PetscFEGetDualSpace(fe, &sp);
1260: PetscDualSpaceGetDimension(sp, &spDim);
1261: for (d = 0; d < spDim; ++d) {
1262: PetscDualSpaceApply(sp, d, geom, numComp, funcs[f], ctx, &values[v]);
1263: v += numComp;
1264: }
1265: }
1266: DMDAVecSetClosure(dm, section, localX, c, values, mode);
1267: }
1268: DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
1269: PetscFree3(v0,J,detJ);
1270: return(0);
1271: }
1275: /*@C
1276: DMDAProjectFunction - This projects the given function into the function space provided.
1278: Input Parameters:
1279: + dm - The DM
1280: . funcs - The coordinate functions to evaluate
1281: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
1282: - mode - The insertion mode for values
1284: Output Parameter:
1285: . X - vector
1287: Level: developer
1289: .seealso: DMDAComputeL2Diff()
1290: @*/
1291: PetscErrorCode DMDAProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
1292: {
1293: Vec localX;
1298: DMGetLocalVector(dm, &localX);
1299: DMDAProjectFunctionLocal(dm, funcs, ctxs, mode, localX);
1300: DMLocalToGlobalBegin(dm, localX, mode, X);
1301: DMLocalToGlobalEnd(dm, localX, mode, X);
1302: DMRestoreLocalVector(dm, &localX);
1303: return(0);
1304: }
1308: /*@C
1309: DMDAComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
1311: Input Parameters:
1312: + dm - The DM
1313: . funcs - The functions to evaluate for each field component
1314: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
1315: - X - The coefficient vector u_h
1317: Output Parameter:
1318: . diff - The diff ||u - u_h||_2
1320: Level: developer
1322: .seealso: DMDAProjectFunction(), DMDAComputeL2GradientDiff()
1323: @*/
1324: PetscErrorCode DMDAComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1325: {
1326: const PetscInt debug = 0;
1327: PetscDS prob;
1328: PetscFE fe;
1329: PetscQuadrature quad;
1330: PetscSection section;
1331: Vec localX;
1332: PetscScalar *funcVal;
1333: PetscReal *coords, *v0, *J, *invJ, detJ;
1334: PetscReal localDiff = 0.0;
1335: PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1336: PetscErrorCode ierr;
1339: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1340: DMGetDS(dm, &prob);
1341: PetscDSGetTotalComponents(prob, &Nc);
1342: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1343: PetscFEGetQuadrature(fe, &quad);
1344: DMGetDefaultSection(dm, §ion);
1345: PetscSectionGetNumFields(section, &numFields);
1346: DMGetLocalVector(dm, &localX);
1347: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1348: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1349: /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1350: PetscMalloc5(Nc,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1351: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1352: for (c = cStart; c < cEnd; ++c) {
1353: PetscScalar *x = NULL;
1354: PetscReal elemDiff = 0.0;
1356: DMDAComputeCellGeometry(dm, c, quad, v0, J, invJ, &detJ);
1357: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1358: DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);
1360: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1361: void * const ctx = ctxs ? ctxs[field] : NULL;
1362: const PetscReal *quadPoints, *quadWeights;
1363: PetscReal *basis;
1364: PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;
1366: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1367: PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1368: PetscFEGetDimension(fe, &numBasisFuncs);
1369: PetscFEGetNumComponents(fe, &numBasisComps);
1370: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
1371: if (debug) {
1372: char title[1024];
1373: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1374: DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
1375: }
1376: for (q = 0; q < numQuadPoints; ++q) {
1377: for (d = 0; d < dim; d++) {
1378: coords[d] = v0[d];
1379: for (e = 0; e < dim; e++) {
1380: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1381: }
1382: }
1383: (*funcs[field])(coords, funcVal, ctx);
1384: for (fc = 0; fc < numBasisComps; ++fc) {
1385: PetscScalar interpolant = 0.0;
1387: for (f = 0; f < numBasisFuncs; ++f) {
1388: const PetscInt fidx = f*numBasisComps+fc;
1389: interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
1390: }
1391: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
1392: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
1393: }
1394: }
1395: comp += numBasisComps;
1396: fieldOffset += numBasisFuncs*numBasisComps;
1397: }
1398: DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1399: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
1400: localDiff += elemDiff;
1401: }
1402: PetscFree5(funcVal,coords,v0,J,invJ);
1403: DMRestoreLocalVector(dm, &localX);
1404: MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
1405: *diff = PetscSqrtReal(*diff);
1406: return(0);
1407: }
1411: /*@C
1412: DMDAComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
1414: Input Parameters:
1415: + dm - The DM
1416: . funcs - The gradient functions to evaluate for each field component
1417: . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null.
1418: . X - The coefficient vector u_h
1419: - n - The vector to project along
1421: Output Parameter:
1422: . diff - The diff ||(grad u - grad u_h) . n||_2
1424: Level: developer
1426: .seealso: DMDAProjectFunction(), DMPlexComputeL2Diff()
1427: @*/
1428: PetscErrorCode DMDAComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
1429: {
1430: const PetscInt debug = 0;
1431: PetscDS prob;
1432: PetscFE fe;
1433: PetscQuadrature quad;
1434: PetscSection section;
1435: Vec localX;
1436: PetscScalar *funcVal, *interpolantVec;
1437: PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
1438: PetscReal localDiff = 0.0;
1439: PetscInt dim, numFields, Nc, cStart, cEnd, c, field, fieldOffset, comp;
1440: PetscErrorCode ierr;
1443: DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);
1444: DMGetDS(dm, &prob);
1445: PetscDSGetTotalComponents(prob, &Nc);
1446: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
1447: PetscFEGetQuadrature(fe, &quad);
1448: DMGetDefaultSection(dm, §ion);
1449: PetscSectionGetNumFields(section, &numFields);
1450: DMGetLocalVector(dm, &localX);
1451: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1452: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1453: /* There are no BC values in DAs right now: DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1454: PetscMalloc7(Nc,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
1455: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
1456: for (c = cStart; c < cEnd; ++c) {
1457: PetscScalar *x = NULL;
1458: PetscReal elemDiff = 0.0;
1460: DMDAComputeCellGeometry(dm, c, quad, v0, J, invJ, &detJ);
1461: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1462: DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);
1464: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
1465: void * const ctx = ctxs ? ctxs[field] : NULL;
1466: const PetscReal *quadPoints, *quadWeights;
1467: PetscReal *basisDer;
1468: PetscInt Nq, Nb, Nc, q, d, e, fc, f, g;
1470: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1471: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
1472: PetscFEGetDimension(fe, &Nb);
1473: PetscFEGetNumComponents(fe, &Nc);
1474: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1475: if (debug) {
1476: char title[1024];
1477: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1478: DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1479: }
1480: for (q = 0; q < Nq; ++q) {
1481: for (d = 0; d < dim; d++) {
1482: coords[d] = v0[d];
1483: for (e = 0; e < dim; e++) {
1484: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1485: }
1486: }
1487: (*funcs[field])(coords, n, funcVal, ctx);
1488: for (fc = 0; fc < Nc; ++fc) {
1489: PetscScalar interpolant = 0.0;
1491: for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
1492: for (f = 0; f < Nb; ++f) {
1493: const PetscInt fidx = f*Nc+fc;
1495: for (d = 0; d < dim; ++d) {
1496: realSpaceDer[d] = 0.0;
1497: for (g = 0; g < dim; ++g) {
1498: realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Nc+fidx)*dim+g];
1499: }
1500: interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
1501: }
1502: }
1503: for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
1504: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
1505: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
1506: }
1507: }
1508: comp += Nc;
1509: fieldOffset += Nb*Nc;
1510: }
1511: DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1512: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
1513: localDiff += elemDiff;
1514: }
1515: PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
1516: DMRestoreLocalVector(dm, &localX);
1517: MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
1518: *diff = PetscSqrtReal(*diff);
1519: return(0);
1520: }
1522: /* ------------------------------------------------------------------- */
1526: /*@C
1527: DMDAGetArray - Gets a work array for a DMDA
1529: Input Parameter:
1530: + da - information about my local patch
1531: - ghosted - do you want arrays for the ghosted or nonghosted patch
1533: Output Parameters:
1534: . vptr - array data structured
1536: Note: The vector values are NOT initialized and may have garbage in them, so you may need
1537: to zero them.
1539: Level: advanced
1541: .seealso: DMDARestoreArray()
1543: @*/
1544: PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
1545: {
1547: PetscInt j,i,xs,ys,xm,ym,zs,zm;
1548: char *iarray_start;
1549: void **iptr = (void**)vptr;
1550: DM_DA *dd = (DM_DA*)da->data;
1554: if (ghosted) {
1555: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1556: if (dd->arrayghostedin[i]) {
1557: *iptr = dd->arrayghostedin[i];
1558: iarray_start = (char*)dd->startghostedin[i];
1559: dd->arrayghostedin[i] = NULL;
1560: dd->startghostedin[i] = NULL;
1562: goto done;
1563: }
1564: }
1565: xs = dd->Xs;
1566: ys = dd->Ys;
1567: zs = dd->Zs;
1568: xm = dd->Xe-dd->Xs;
1569: ym = dd->Ye-dd->Ys;
1570: zm = dd->Ze-dd->Zs;
1571: } else {
1572: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1573: if (dd->arrayin[i]) {
1574: *iptr = dd->arrayin[i];
1575: iarray_start = (char*)dd->startin[i];
1576: dd->arrayin[i] = NULL;
1577: dd->startin[i] = NULL;
1579: goto done;
1580: }
1581: }
1582: xs = dd->xs;
1583: ys = dd->ys;
1584: zs = dd->zs;
1585: xm = dd->xe-dd->xs;
1586: ym = dd->ye-dd->ys;
1587: zm = dd->ze-dd->zs;
1588: }
1590: switch (dd->dim) {
1591: case 1: {
1592: void *ptr;
1594: PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);
1596: ptr = (void*)(iarray_start - xs*sizeof(PetscScalar));
1597: *iptr = (void*)ptr;
1598: break;
1599: }
1600: case 2: {
1601: void **ptr;
1603: PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);
1605: ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
1606: for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
1607: *iptr = (void*)ptr;
1608: break;
1609: }
1610: case 3: {
1611: void ***ptr,**bptr;
1613: PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);
1615: ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
1616: bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
1617: for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
1618: for (i=zs; i<zs+zm; i++) {
1619: for (j=ys; j<ys+ym; j++) {
1620: ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
1621: }
1622: }
1623: *iptr = (void*)ptr;
1624: break;
1625: }
1626: default:
1627: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
1628: }
1630: done:
1631: /* add arrays to the checked out list */
1632: if (ghosted) {
1633: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1634: if (!dd->arrayghostedout[i]) {
1635: dd->arrayghostedout[i] = *iptr;
1636: dd->startghostedout[i] = iarray_start;
1637: break;
1638: }
1639: }
1640: } else {
1641: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1642: if (!dd->arrayout[i]) {
1643: dd->arrayout[i] = *iptr;
1644: dd->startout[i] = iarray_start;
1645: break;
1646: }
1647: }
1648: }
1649: return(0);
1650: }
1654: /*@C
1655: DMDARestoreArray - Restores an array of derivative types for a DMDA
1657: Input Parameter:
1658: + da - information about my local patch
1659: . ghosted - do you want arrays for the ghosted or nonghosted patch
1660: - vptr - array data structured to be passed to ad_FormFunctionLocal()
1662: Level: advanced
1664: .seealso: DMDAGetArray()
1666: @*/
1667: PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
1668: {
1669: PetscInt i;
1670: void **iptr = (void**)vptr,*iarray_start = 0;
1671: DM_DA *dd = (DM_DA*)da->data;
1675: if (ghosted) {
1676: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1677: if (dd->arrayghostedout[i] == *iptr) {
1678: iarray_start = dd->startghostedout[i];
1679: dd->arrayghostedout[i] = NULL;
1680: dd->startghostedout[i] = NULL;
1681: break;
1682: }
1683: }
1684: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1685: if (!dd->arrayghostedin[i]) {
1686: dd->arrayghostedin[i] = *iptr;
1687: dd->startghostedin[i] = iarray_start;
1688: break;
1689: }
1690: }
1691: } else {
1692: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1693: if (dd->arrayout[i] == *iptr) {
1694: iarray_start = dd->startout[i];
1695: dd->arrayout[i] = NULL;
1696: dd->startout[i] = NULL;
1697: break;
1698: }
1699: }
1700: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
1701: if (!dd->arrayin[i]) {
1702: dd->arrayin[i] = *iptr;
1703: dd->startin[i] = iarray_start;
1704: break;
1705: }
1706: }
1707: }
1708: return(0);
1709: }