Actual source code: dalocal.c
petsc-3.4.5 2014-06-29
2: /*
3: Code for manipulating distributed regular arrays in parallel.
4: */
6: #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/
7: #include <petscsf.h>
9: /*
10: This allows the DMDA vectors to properly tell MATLAB their dimensions
11: */
12: #if defined(PETSC_HAVE_MATLAB_ENGINE)
13: #include <engine.h> /* MATLAB include file */
14: #include <mex.h> /* MATLAB include file */
17: static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine)
18: {
20: PetscInt n,m;
21: Vec vec = (Vec)obj;
22: PetscScalar *array;
23: mxArray *mat;
24: DM da;
27: VecGetDM(vec, &da);
28: if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA");
29: DMDAGetGhostCorners(da,0,0,0,&m,&n,0);
31: VecGetArray(vec,&array);
32: #if !defined(PETSC_USE_COMPLEX)
33: mat = mxCreateDoubleMatrix(m,n,mxREAL);
34: #else
35: mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
36: #endif
37: PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));
38: PetscObjectName(obj);
39: engPutVariable((Engine*)mengine,obj->name,mat);
41: VecRestoreArray(vec,&array);
42: return(0);
43: }
44: #endif
49: PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g)
50: {
52: DM_DA *dd = (DM_DA*)da->data;
57: if (da->defaultSection) {
58: DMCreateLocalVector_Section_Private(da,g);
59: } else {
60: VecCreate(PETSC_COMM_SELF,g);
61: VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);
62: VecSetBlockSize(*g,dd->w);
63: VecSetType(*g,da->vectype);
64: VecSetDM(*g, da);
65: #if defined(PETSC_HAVE_MATLAB_ENGINE)
66: if (dd->w == 1 && dd->dim == 2) {
67: PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);
68: }
69: #endif
70: }
71: return(0);
72: }
76: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCells)
77: {
78: DM_DA *da = (DM_DA*) dm->data;
79: const PetscInt dim = da->dim;
80: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
81: const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
84: if (numCells) {
86: *numCells = nC;
87: }
88: return(0);
89: }
93: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
94: {
95: DM_DA *da = (DM_DA*) dm->data;
96: const PetscInt dim = da->dim;
97: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
98: const PetscInt nVx = mx+1;
99: const PetscInt nVy = dim > 1 ? (my+1) : 1;
100: const PetscInt nVz = dim > 2 ? (mz+1) : 1;
101: const PetscInt nV = nVx*nVy*nVz;
104: if (numVerticesX) {
106: *numVerticesX = nVx;
107: }
108: if (numVerticesY) {
110: *numVerticesY = nVy;
111: }
112: if (numVerticesZ) {
114: *numVerticesZ = nVz;
115: }
116: if (numVertices) {
118: *numVertices = nV;
119: }
120: return(0);
121: }
125: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
126: {
127: DM_DA *da = (DM_DA*) dm->data;
128: const PetscInt dim = da->dim;
129: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
130: const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
131: const PetscInt nXF = (mx+1)*nxF;
132: const PetscInt nyF = mx*(dim > 2 ? mz : 1);
133: const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0;
134: const PetscInt nzF = mx*(dim > 1 ? my : 0);
135: const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0;
138: if (numXFacesX) {
140: *numXFacesX = nxF;
141: }
142: if (numXFaces) {
144: *numXFaces = nXF;
145: }
146: if (numYFacesY) {
148: *numYFacesY = nyF;
149: }
150: if (numYFaces) {
152: *numYFaces = nYF;
153: }
154: if (numZFacesZ) {
156: *numZFacesZ = nzF;
157: }
158: if (numZFaces) {
160: *numZFaces = nZF;
161: }
162: return(0);
163: }
167: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
168: {
169: DM_DA *da = (DM_DA*) dm->data;
170: const PetscInt dim = da->dim;
171: PetscInt nC, nV, nXF, nYF, nZF;
177: DMDAGetNumCells(dm, &nC);
178: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
179: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
180: if (height == 0) {
181: /* Cells */
182: if (pStart) *pStart = 0;
183: if (pEnd) *pEnd = nC;
184: } else if (height == 1) {
185: /* Faces */
186: if (pStart) *pStart = nC+nV;
187: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
188: } else if (height == dim) {
189: /* Vertices */
190: if (pStart) *pStart = nC;
191: if (pEnd) *pEnd = nC+nV;
192: } else if (height < 0) {
193: /* All points */
194: if (pStart) *pStart = 0;
195: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
196: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
197: return(0);
198: }
202: /*@C
203: DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with
204: different numbers of dofs on vertices, cells, and faces in each direction.
206: Input Parameters:
207: + dm- The DMDA
208: . numFields - The number of fields
209: . numComp - The number of components in each field, or NULL for 1
210: . numVertexDof - The number of dofs per vertex for each field, or NULL
211: . numFaceDof - The number of dofs per face for each field and direction, or NULL
212: - numCellDof - The number of dofs per cell for each field, or NULL
214: Level: developer
216: Note:
217: The default DMDA numbering is as follows:
219: - Cells: [0, nC)
220: - Vertices: [nC, nC+nV)
221: - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir
222: - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir
223: - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir
225: We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment.
226: @*/
227: PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numVertexDof[], PetscInt numFaceDof[], PetscInt numCellDof[])
228: {
229: DM_DA *da = (DM_DA*) dm->data;
230: const PetscInt dim = da->dim;
231: PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0};
232: PetscSF sf;
233: PetscMPIInt rank;
234: const PetscMPIInt *neighbors;
235: PetscInt *localPoints;
236: PetscSFNode *remotePoints;
237: PetscInt nleaves = 0, nleavesCheck = 0, nL = 0;
238: PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF;
239: PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd;
240: PetscInt f, v, c, xf, yf, zf, xn, yn, zn;
241: PetscErrorCode ierr;
245: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
246: DMDAGetNumCells(dm, &nC);
247: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
248: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
249: DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);
250: DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);
251: DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);
252: DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);
253: xfStart = vEnd; xfEnd = xfStart+nXF;
254: yfStart = xfEnd; yfEnd = yfStart+nYF;
255: zfStart = yfEnd; zfEnd = zfStart+nZF;
256: if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd);
257: /* Create local section */
258: DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);
259: for (f = 0; f < numFields; ++f) {
260: if (numVertexDof) numVertexTotDof += numVertexDof[f];
261: if (numCellDof) numCellTotDof += numCellDof[f];
262: if (numFaceDof) {
263: numFaceTotDof[0] += numFaceDof[f*dim+0];
264: numFaceTotDof[1] += dim > 1 ? numFaceDof[f*dim+1] : 0;
265: numFaceTotDof[2] += dim > 2 ? numFaceDof[f*dim+2] : 0;
266: }
267: }
268: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &dm->defaultSection);
269: if (numFields > 1) {
270: PetscSectionSetNumFields(dm->defaultSection, numFields);
271: for (f = 0; f < numFields; ++f) {
272: const char *name;
274: DMDAGetFieldName(dm, f, &name);
275: PetscSectionSetFieldName(dm->defaultSection, f, name);
276: if (numComp) {
277: PetscSectionSetFieldComponents(dm->defaultSection, f, numComp[f]);
278: }
279: }
280: } else {
281: numFields = 0;
282: }
283: PetscSectionSetChart(dm->defaultSection, pStart, pEnd);
284: if (numVertexDof) {
285: for (v = vStart; v < vEnd; ++v) {
286: for (f = 0; f < numFields; ++f) {
287: PetscSectionSetFieldDof(dm->defaultSection, v, f, numVertexDof[f]);
288: }
289: PetscSectionSetDof(dm->defaultSection, v, numVertexTotDof);
290: }
291: }
292: if (numFaceDof) {
293: for (xf = xfStart; xf < xfEnd; ++xf) {
294: for (f = 0; f < numFields; ++f) {
295: PetscSectionSetFieldDof(dm->defaultSection, xf, f, numFaceDof[f*dim+0]);
296: }
297: PetscSectionSetDof(dm->defaultSection, xf, numFaceTotDof[0]);
298: }
299: for (yf = yfStart; yf < yfEnd; ++yf) {
300: for (f = 0; f < numFields; ++f) {
301: PetscSectionSetFieldDof(dm->defaultSection, yf, f, numFaceDof[f*dim+1]);
302: }
303: PetscSectionSetDof(dm->defaultSection, yf, numFaceTotDof[1]);
304: }
305: for (zf = zfStart; zf < zfEnd; ++zf) {
306: for (f = 0; f < numFields; ++f) {
307: PetscSectionSetFieldDof(dm->defaultSection, zf, f, numFaceDof[f*dim+2]);
308: }
309: PetscSectionSetDof(dm->defaultSection, zf, numFaceTotDof[2]);
310: }
311: }
312: if (numCellDof) {
313: for (c = cStart; c < cEnd; ++c) {
314: for (f = 0; f < numFields; ++f) {
315: PetscSectionSetFieldDof(dm->defaultSection, c, f, numCellDof[f]);
316: }
317: PetscSectionSetDof(dm->defaultSection, c, numCellTotDof);
318: }
319: }
320: PetscSectionSetUp(dm->defaultSection);
321: /* Create mesh point SF */
322: DMDAGetNeighbors(dm, &neighbors);
323: for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
324: for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
325: for (xn = 0; xn < 3; ++xn) {
326: const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
327: const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
329: if (neighbor >= 0 && neighbor < rank) {
330: nleaves += (!xp ? nVx : 1) * (!yp ? nVy : 1) * (!zp ? nVz : 1); /* vertices */
331: if (xp && !yp && !zp) {
332: nleaves += nxF; /* x faces */
333: } else if (yp && !zp && !xp) {
334: nleaves += nyF; /* y faces */
335: } else if (zp && !xp && !yp) {
336: nleaves += nzF; /* z faces */
337: }
338: }
339: }
340: }
341: }
342: PetscMalloc2(nleaves,PetscInt,&localPoints,nleaves,PetscSFNode,&remotePoints);
343: for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) {
344: for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) {
345: for (xn = 0; xn < 3; ++xn) {
346: const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0;
347: const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn];
348: PetscInt xv, yv, zv;
350: if (neighbor >= 0 && neighbor < rank) {
351: if (xp < 0) { /* left */
352: if (yp < 0) { /* bottom */
353: if (zp < 0) { /* back */
354: const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC;
355: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
356: nleavesCheck += 1; /* left bottom back vertex */
358: localPoints[nL] = localVertex;
359: remotePoints[nL].rank = neighbor;
360: remotePoints[nL].index = remoteVertex;
361: ++nL;
362: } else if (zp > 0) { /* front */
363: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC;
364: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
365: nleavesCheck += 1; /* left bottom front vertex */
367: localPoints[nL] = localVertex;
368: remotePoints[nL].rank = neighbor;
369: remotePoints[nL].index = remoteVertex;
370: ++nL;
371: } else {
372: nleavesCheck += nVz; /* left bottom vertices */
373: for (zv = 0; zv < nVz; ++zv, ++nL) {
374: const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC;
375: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
377: localPoints[nL] = localVertex;
378: remotePoints[nL].rank = neighbor;
379: remotePoints[nL].index = remoteVertex;
380: }
381: }
382: } else if (yp > 0) { /* top */
383: if (zp < 0) { /* back */
384: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC;
385: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
386: nleavesCheck += 1; /* left top back vertex */
388: localPoints[nL] = localVertex;
389: remotePoints[nL].rank = neighbor;
390: remotePoints[nL].index = remoteVertex;
391: ++nL;
392: } else if (zp > 0) { /* front */
393: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC;
394: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
395: nleavesCheck += 1; /* left top front vertex */
397: localPoints[nL] = localVertex;
398: remotePoints[nL].rank = neighbor;
399: remotePoints[nL].index = remoteVertex;
400: ++nL;
401: } else {
402: nleavesCheck += nVz; /* left top vertices */
403: for (zv = 0; zv < nVz; ++zv, ++nL) {
404: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC;
405: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
407: localPoints[nL] = localVertex;
408: remotePoints[nL].rank = neighbor;
409: remotePoints[nL].index = remoteVertex;
410: }
411: }
412: } else {
413: if (zp < 0) { /* back */
414: nleavesCheck += nVy; /* left back vertices */
415: for (yv = 0; yv < nVy; ++yv, ++nL) {
416: const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC;
417: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
419: localPoints[nL] = localVertex;
420: remotePoints[nL].rank = neighbor;
421: remotePoints[nL].index = remoteVertex;
422: }
423: } else if (zp > 0) { /* front */
424: nleavesCheck += nVy; /* left front vertices */
425: for (yv = 0; yv < nVy; ++yv, ++nL) {
426: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC;
427: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
429: localPoints[nL] = localVertex;
430: remotePoints[nL].rank = neighbor;
431: remotePoints[nL].index = remoteVertex;
432: }
433: } else {
434: nleavesCheck += nVy*nVz; /* left vertices */
435: for (zv = 0; zv < nVz; ++zv) {
436: for (yv = 0; yv < nVy; ++yv, ++nL) {
437: const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC;
438: const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */
440: localPoints[nL] = localVertex;
441: remotePoints[nL].rank = neighbor;
442: remotePoints[nL].index = remoteVertex;
443: }
444: }
445: nleavesCheck += nxF; /* left faces */
446: for (xf = 0; xf < nxF; ++xf, ++nL) {
447: /* THIS IS WRONG */
448: const PetscInt localFace = 0 + nC+nV;
449: const PetscInt remoteFace = 0 + nC+nV;
451: localPoints[nL] = localFace;
452: remotePoints[nL].rank = neighbor;
453: remotePoints[nL].index = remoteFace;
454: }
455: }
456: }
457: } else if (xp > 0) { /* right */
458: if (yp < 0) { /* bottom */
459: if (zp < 0) { /* back */
460: const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC;
461: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
462: nleavesCheck += 1; /* right bottom back vertex */
464: localPoints[nL] = localVertex;
465: remotePoints[nL].rank = neighbor;
466: remotePoints[nL].index = remoteVertex;
467: ++nL;
468: } else if (zp > 0) { /* front */
469: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC;
470: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
471: nleavesCheck += 1; /* right bottom front vertex */
473: localPoints[nL] = localVertex;
474: remotePoints[nL].rank = neighbor;
475: remotePoints[nL].index = remoteVertex;
476: ++nL;
477: } else {
478: nleavesCheck += nVz; /* right bottom vertices */
479: for (zv = 0; zv < nVz; ++zv, ++nL) {
480: const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC;
481: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
483: localPoints[nL] = localVertex;
484: remotePoints[nL].rank = neighbor;
485: remotePoints[nL].index = remoteVertex;
486: }
487: }
488: } else if (yp > 0) { /* top */
489: if (zp < 0) { /* back */
490: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC;
491: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
492: nleavesCheck += 1; /* right top back vertex */
494: localPoints[nL] = localVertex;
495: remotePoints[nL].rank = neighbor;
496: remotePoints[nL].index = remoteVertex;
497: ++nL;
498: } else if (zp > 0) { /* front */
499: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC;
500: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
501: nleavesCheck += 1; /* right top front vertex */
503: localPoints[nL] = localVertex;
504: remotePoints[nL].rank = neighbor;
505: remotePoints[nL].index = remoteVertex;
506: ++nL;
507: } else {
508: nleavesCheck += nVz; /* right top vertices */
509: for (zv = 0; zv < nVz; ++zv, ++nL) {
510: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC;
511: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
513: localPoints[nL] = localVertex;
514: remotePoints[nL].rank = neighbor;
515: remotePoints[nL].index = remoteVertex;
516: }
517: }
518: } else {
519: if (zp < 0) { /* back */
520: nleavesCheck += nVy; /* right back vertices */
521: for (yv = 0; yv < nVy; ++yv, ++nL) {
522: const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC;
523: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
525: localPoints[nL] = localVertex;
526: remotePoints[nL].rank = neighbor;
527: remotePoints[nL].index = remoteVertex;
528: }
529: } else if (zp > 0) { /* front */
530: nleavesCheck += nVy; /* right front vertices */
531: for (yv = 0; yv < nVy; ++yv, ++nL) {
532: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC;
533: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
535: localPoints[nL] = localVertex;
536: remotePoints[nL].rank = neighbor;
537: remotePoints[nL].index = remoteVertex;
538: }
539: } else {
540: nleavesCheck += nVy*nVz; /* right vertices */
541: for (zv = 0; zv < nVz; ++zv) {
542: for (yv = 0; yv < nVy; ++yv, ++nL) {
543: const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC;
544: const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */
546: localPoints[nL] = localVertex;
547: remotePoints[nL].rank = neighbor;
548: remotePoints[nL].index = remoteVertex;
549: }
550: }
551: nleavesCheck += nxF; /* right faces */
552: for (xf = 0; xf < nxF; ++xf, ++nL) {
553: /* THIS IS WRONG */
554: const PetscInt localFace = 0 + nC+nV;
555: const PetscInt remoteFace = 0 + nC+nV;
557: localPoints[nL] = localFace;
558: remotePoints[nL].rank = neighbor;
559: remotePoints[nL].index = remoteFace;
560: }
561: }
562: }
563: } else {
564: if (yp < 0) { /* bottom */
565: if (zp < 0) { /* back */
566: nleavesCheck += nVx; /* bottom back vertices */
567: for (xv = 0; xv < nVx; ++xv, ++nL) {
568: const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC;
569: const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
571: localPoints[nL] = localVertex;
572: remotePoints[nL].rank = neighbor;
573: remotePoints[nL].index = remoteVertex;
574: }
575: } else if (zp > 0) { /* front */
576: nleavesCheck += nVx; /* bottom front vertices */
577: for (xv = 0; xv < nVx; ++xv, ++nL) {
578: const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC;
579: const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
581: localPoints[nL] = localVertex;
582: remotePoints[nL].rank = neighbor;
583: remotePoints[nL].index = remoteVertex;
584: }
585: } else {
586: nleavesCheck += nVx*nVz; /* bottom vertices */
587: for (zv = 0; zv < nVz; ++zv) {
588: for (xv = 0; xv < nVx; ++xv, ++nL) {
589: const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC;
590: const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
592: localPoints[nL] = localVertex;
593: remotePoints[nL].rank = neighbor;
594: remotePoints[nL].index = remoteVertex;
595: }
596: }
597: nleavesCheck += nyF; /* bottom faces */
598: for (yf = 0; yf < nyF; ++yf, ++nL) {
599: /* THIS IS WRONG */
600: const PetscInt localFace = 0 + nC+nV;
601: const PetscInt remoteFace = 0 + nC+nV;
603: localPoints[nL] = localFace;
604: remotePoints[nL].rank = neighbor;
605: remotePoints[nL].index = remoteFace;
606: }
607: }
608: } else if (yp > 0) { /* top */
609: if (zp < 0) { /* back */
610: nleavesCheck += nVx; /* top back vertices */
611: for (xv = 0; xv < nVx; ++xv, ++nL) {
612: const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC;
613: const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
615: localPoints[nL] = localVertex;
616: remotePoints[nL].rank = neighbor;
617: remotePoints[nL].index = remoteVertex;
618: }
619: } else if (zp > 0) { /* front */
620: nleavesCheck += nVx; /* top front vertices */
621: for (xv = 0; xv < nVx; ++xv, ++nL) {
622: const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC;
623: const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
625: localPoints[nL] = localVertex;
626: remotePoints[nL].rank = neighbor;
627: remotePoints[nL].index = remoteVertex;
628: }
629: } else {
630: nleavesCheck += nVx*nVz; /* top vertices */
631: for (zv = 0; zv < nVz; ++zv) {
632: for (xv = 0; xv < nVx; ++xv, ++nL) {
633: const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC;
634: const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
636: localPoints[nL] = localVertex;
637: remotePoints[nL].rank = neighbor;
638: remotePoints[nL].index = remoteVertex;
639: }
640: }
641: nleavesCheck += nyF; /* top faces */
642: for (yf = 0; yf < nyF; ++yf, ++nL) {
643: /* THIS IS WRONG */
644: const PetscInt localFace = 0 + nC+nV;
645: const PetscInt remoteFace = 0 + nC+nV;
647: localPoints[nL] = localFace;
648: remotePoints[nL].rank = neighbor;
649: remotePoints[nL].index = remoteFace;
650: }
651: }
652: } else {
653: if (zp < 0) { /* back */
654: nleavesCheck += nVx*nVy; /* back vertices */
655: for (yv = 0; yv < nVy; ++yv) {
656: for (xv = 0; xv < nVx; ++xv, ++nL) {
657: const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC;
658: const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
660: localPoints[nL] = localVertex;
661: remotePoints[nL].rank = neighbor;
662: remotePoints[nL].index = remoteVertex;
663: }
664: }
665: nleavesCheck += nzF; /* back faces */
666: for (zf = 0; zf < nzF; ++zf, ++nL) {
667: /* THIS IS WRONG */
668: const PetscInt localFace = 0 + nC+nV;
669: const PetscInt remoteFace = 0 + nC+nV;
671: localPoints[nL] = localFace;
672: remotePoints[nL].rank = neighbor;
673: remotePoints[nL].index = remoteFace;
674: }
675: } else if (zp > 0) { /* front */
676: nleavesCheck += nVx*nVy; /* front vertices */
677: for (yv = 0; yv < nVy; ++yv) {
678: for (xv = 0; xv < nVx; ++xv, ++nL) {
679: const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC;
680: const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */
682: localPoints[nL] = localVertex;
683: remotePoints[nL].rank = neighbor;
684: remotePoints[nL].index = remoteVertex;
685: }
686: }
687: nleavesCheck += nzF; /* front faces */
688: for (zf = 0; zf < nzF; ++zf, ++nL) {
689: /* THIS IS WRONG */
690: const PetscInt localFace = 0 + nC+nV;
691: const PetscInt remoteFace = 0 + nC+nV;
693: localPoints[nL] = localFace;
694: remotePoints[nL].rank = neighbor;
695: remotePoints[nL].index = remoteFace;
696: }
697: } else {
698: /* Nothing is shared from the interior */
699: }
700: }
701: }
702: }
703: }
704: }
705: }
706: /* TODO: Remove duplication in leaf determination */
707: if (nleaves != nleavesCheck) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck);
708: PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);
709: PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
710: /* Create global section */
711: PetscSectionCreateGlobalSection(dm->defaultSection, sf, PETSC_FALSE, &dm->defaultGlobalSection);
712: PetscSFDestroy(&sf);
713: /* Create default SF */
714: DMCreateDefaultSF(dm, dm->defaultSection, dm->defaultGlobalSection);
715: return(0);
716: }
718: /* ------------------------------------------------------------------- */
722: /*@C
723: DMDAGetArray - Gets a work array for a DMDA
725: Input Parameter:
726: + da - information about my local patch
727: - ghosted - do you want arrays for the ghosted or nonghosted patch
729: Output Parameters:
730: . vptr - array data structured
732: Note: The vector values are NOT initialized and may have garbage in them, so you may need
733: to zero them.
735: Level: advanced
737: .seealso: DMDARestoreArray()
739: @*/
740: PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
741: {
743: PetscInt j,i,xs,ys,xm,ym,zs,zm;
744: char *iarray_start;
745: void **iptr = (void**)vptr;
746: DM_DA *dd = (DM_DA*)da->data;
750: if (ghosted) {
751: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
752: if (dd->arrayghostedin[i]) {
753: *iptr = dd->arrayghostedin[i];
754: iarray_start = (char*)dd->startghostedin[i];
755: dd->arrayghostedin[i] = NULL;
756: dd->startghostedin[i] = NULL;
758: goto done;
759: }
760: }
761: xs = dd->Xs;
762: ys = dd->Ys;
763: zs = dd->Zs;
764: xm = dd->Xe-dd->Xs;
765: ym = dd->Ye-dd->Ys;
766: zm = dd->Ze-dd->Zs;
767: } else {
768: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
769: if (dd->arrayin[i]) {
770: *iptr = dd->arrayin[i];
771: iarray_start = (char*)dd->startin[i];
772: dd->arrayin[i] = NULL;
773: dd->startin[i] = NULL;
775: goto done;
776: }
777: }
778: xs = dd->xs;
779: ys = dd->ys;
780: zs = dd->zs;
781: xm = dd->xe-dd->xs;
782: ym = dd->ye-dd->ys;
783: zm = dd->ze-dd->zs;
784: }
786: switch (dd->dim) {
787: case 1: {
788: void *ptr;
790: PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);
792: ptr = (void*)(iarray_start - xs*sizeof(PetscScalar));
793: *iptr = (void*)ptr;
794: break;
795: }
796: case 2: {
797: void **ptr;
799: PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);
801: ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
802: for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
803: *iptr = (void*)ptr;
804: break;
805: }
806: case 3: {
807: void ***ptr,**bptr;
809: PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);
811: ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
812: bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
813: for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
814: for (i=zs; i<zs+zm; i++) {
815: for (j=ys; j<ys+ym; j++) {
816: ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
817: }
818: }
819: *iptr = (void*)ptr;
820: break;
821: }
822: default:
823: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim);
824: }
826: done:
827: /* add arrays to the checked out list */
828: if (ghosted) {
829: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
830: if (!dd->arrayghostedout[i]) {
831: dd->arrayghostedout[i] = *iptr;
832: dd->startghostedout[i] = iarray_start;
833: break;
834: }
835: }
836: } else {
837: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
838: if (!dd->arrayout[i]) {
839: dd->arrayout[i] = *iptr;
840: dd->startout[i] = iarray_start;
841: break;
842: }
843: }
844: }
845: return(0);
846: }
850: /*@C
851: DMDARestoreArray - Restores an array of derivative types for a DMDA
853: Input Parameter:
854: + da - information about my local patch
855: . ghosted - do you want arrays for the ghosted or nonghosted patch
856: - vptr - array data structured to be passed to ad_FormFunctionLocal()
858: Level: advanced
860: .seealso: DMDAGetArray()
862: @*/
863: PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
864: {
865: PetscInt i;
866: void **iptr = (void**)vptr,*iarray_start = 0;
867: DM_DA *dd = (DM_DA*)da->data;
871: if (ghosted) {
872: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
873: if (dd->arrayghostedout[i] == *iptr) {
874: iarray_start = dd->startghostedout[i];
875: dd->arrayghostedout[i] = NULL;
876: dd->startghostedout[i] = NULL;
877: break;
878: }
879: }
880: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
881: if (!dd->arrayghostedin[i]) {
882: dd->arrayghostedin[i] = *iptr;
883: dd->startghostedin[i] = iarray_start;
884: break;
885: }
886: }
887: } else {
888: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
889: if (dd->arrayout[i] == *iptr) {
890: iarray_start = dd->startout[i];
891: dd->arrayout[i] = NULL;
892: dd->startout[i] = NULL;
893: break;
894: }
895: }
896: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
897: if (!dd->arrayin[i]) {
898: dd->arrayin[i] = *iptr;
899: dd->startin[i] = iarray_start;
900: break;
901: }
902: }
903: }
904: return(0);
905: }