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