Actual source code: dalocal.c
2: /*
3: Code for manipulating distributed regular arrays in parallel.
4: */
6: #include <petsc/private/dmdaimpl.h>
7: #include <petscbt.h>
8: #include <petscsf.h>
9: #include <petscds.h>
10: #include <petscfe.h>
12: /*
13: This allows the DMDA vectors to properly tell MATLAB their dimensions
14: */
15: #if defined(PETSC_HAVE_MATLAB)
16: #include <engine.h> /* MATLAB include file */
17: #include <mex.h> /* MATLAB include file */
18: static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj, void *mengine)
19: {
20: PetscInt n, m;
21: Vec vec = (Vec)obj;
22: PetscScalar *array;
23: mxArray *mat;
24: DM da;
26: VecGetDM(vec, &da);
28: DMDAGetGhostCorners(da, 0, 0, 0, &m, &n, 0);
30: VecGetArray(vec, &array);
31: #if !defined(PETSC_USE_COMPLEX)
32: mat = mxCreateDoubleMatrix(m, n, mxREAL);
33: #else
34: mat = mxCreateDoubleMatrix(m, n, mxCOMPLEX);
35: #endif
36: PetscArraycpy(mxGetPr(mat), array, n * m);
37: PetscObjectName(obj);
38: engPutVariable((Engine *)mengine, obj->name, mat);
40: VecRestoreArray(vec, &array);
41: return 0;
42: }
43: #endif
45: PetscErrorCode DMCreateLocalVector_DA(DM da, Vec *g)
46: {
47: DM_DA *dd = (DM_DA *)da->data;
51: VecCreate(PETSC_COMM_SELF, g);
52: VecSetSizes(*g, dd->nlocal, PETSC_DETERMINE);
53: VecSetBlockSize(*g, dd->w);
54: VecSetType(*g, da->vectype);
55: if (dd->nlocal < da->bind_below) {
56: VecSetBindingPropagates(*g, PETSC_TRUE);
57: VecBindToCPU(*g, PETSC_TRUE);
58: }
59: VecSetDM(*g, da);
60: #if defined(PETSC_HAVE_MATLAB)
61: if (dd->w == 1 && da->dim == 2) PetscObjectComposeFunction((PetscObject)*g, "PetscMatlabEnginePut_C", VecMatlabEnginePut_DA2d);
62: #endif
63: return 0;
64: }
66: /*@
67: DMDAGetNumCells - Get the number of cells in the local piece of the `DMDA`. This includes ghost cells.
69: Input Parameter:
70: . dm - The `DMDA` object
72: Output Parameters:
73: + numCellsX - The number of local cells in the x-direction
74: . numCellsY - The number of local cells in the y-direction
75: . numCellsZ - The number of local cells in the z-direction
76: - numCells - The number of local cells
78: Level: developer
80: .seealso: `DM`, `DMDA`, `DMDAGetCellPoint()`
81: @*/
82: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
83: {
84: DM_DA *da = (DM_DA *)dm->data;
85: const PetscInt dim = dm->dim;
86: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
87: const PetscInt nC = (mx) * (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
90: if (numCellsX) {
92: *numCellsX = mx;
93: }
94: if (numCellsY) {
96: *numCellsY = my;
97: }
98: if (numCellsZ) {
100: *numCellsZ = mz;
101: }
102: if (numCells) {
104: *numCells = nC;
105: }
106: return 0;
107: }
109: /*@
110: DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the `DMDA`
112: Input Parameters:
113: + dm - The `DMDA` object
114: - i,j,k - The global indices for the cell
116: Output Parameters:
117: . point - The local `DM` point
119: Level: developer
121: .seealso: `DM`, `DMDA`, `DMDAGetNumCells()`
122: @*/
123: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
124: {
125: const PetscInt dim = dm->dim;
126: DMDALocalInfo info;
130: DMDAGetLocalInfo(dm, &info);
134: *point = i + (dim > 1 ? (j + (dim > 2 ? k * info.gym : 0)) * info.gxm : 0);
135: return 0;
136: }
138: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
139: {
140: DM_DA *da = (DM_DA *)dm->data;
141: const PetscInt dim = dm->dim;
142: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
143: const PetscInt nVx = mx + 1;
144: const PetscInt nVy = dim > 1 ? (my + 1) : 1;
145: const PetscInt nVz = dim > 2 ? (mz + 1) : 1;
146: const PetscInt nV = nVx * nVy * nVz;
148: if (numVerticesX) {
150: *numVerticesX = nVx;
151: }
152: if (numVerticesY) {
154: *numVerticesY = nVy;
155: }
156: if (numVerticesZ) {
158: *numVerticesZ = nVz;
159: }
160: if (numVertices) {
162: *numVertices = nV;
163: }
164: return 0;
165: }
167: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
168: {
169: DM_DA *da = (DM_DA *)dm->data;
170: const PetscInt dim = dm->dim;
171: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
172: const PetscInt nxF = (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
173: const PetscInt nXF = (mx + 1) * nxF;
174: const PetscInt nyF = mx * (dim > 2 ? mz : 1);
175: const PetscInt nYF = dim > 1 ? (my + 1) * nyF : 0;
176: const PetscInt nzF = mx * (dim > 1 ? my : 0);
177: const PetscInt nZF = dim > 2 ? (mz + 1) * nzF : 0;
179: if (numXFacesX) {
181: *numXFacesX = nxF;
182: }
183: if (numXFaces) {
185: *numXFaces = nXF;
186: }
187: if (numYFacesY) {
189: *numYFacesY = nyF;
190: }
191: if (numYFaces) {
193: *numYFaces = nYF;
194: }
195: if (numZFacesZ) {
197: *numZFacesZ = nzF;
198: }
199: if (numZFaces) {
201: *numZFaces = nZF;
202: }
203: return 0;
204: }
206: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
207: {
208: const PetscInt dim = dm->dim;
209: PetscInt nC, nV, nXF, nYF, nZF;
213: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
214: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
215: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
216: if (height == 0) {
217: /* Cells */
218: if (pStart) *pStart = 0;
219: if (pEnd) *pEnd = nC;
220: } else if (height == 1) {
221: /* Faces */
222: if (pStart) *pStart = nC + nV;
223: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
224: } else if (height == dim) {
225: /* Vertices */
226: if (pStart) *pStart = nC;
227: if (pEnd) *pEnd = nC + nV;
228: } else if (height < 0) {
229: /* All points */
230: if (pStart) *pStart = 0;
231: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
232: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %" PetscInt_FMT " in the DA", height);
233: return 0;
234: }
236: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
237: {
238: const PetscInt dim = dm->dim;
239: 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 (depth == dim) {
247: /* Cells */
248: if (pStart) *pStart = 0;
249: if (pEnd) *pEnd = nC;
250: } else if (depth == dim - 1) {
251: /* Faces */
252: if (pStart) *pStart = nC + nV;
253: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
254: } else if (depth == 0) {
255: /* Vertices */
256: if (pStart) *pStart = nC;
257: if (pEnd) *pEnd = nC + nV;
258: } else if (depth < 0) {
259: /* All points */
260: if (pStart) *pStart = 0;
261: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
262: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %" PetscInt_FMT " in the DA", depth);
263: return 0;
264: }
266: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
267: {
268: const PetscInt dim = dm->dim;
269: PetscInt nC, nV, nXF, nYF, nZF;
271: *coneSize = 0;
272: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
273: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
274: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
275: switch (dim) {
276: case 2:
277: if (p >= 0) {
278: if (p < nC) {
279: *coneSize = 4;
280: } else if (p < nC + nV) {
281: *coneSize = 0;
282: } else if (p < nC + nV + nXF + nYF + nZF) {
283: *coneSize = 2;
284: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", p, nC + nV + nXF + nYF + nZF);
285: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %" PetscInt_FMT " is invalid", p);
286: break;
287: case 3:
288: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
289: }
290: return 0;
291: }
293: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
294: {
295: const PetscInt dim = dm->dim;
296: PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
298: if (!cone) DMGetWorkArray(dm, 6, MPIU_INT, cone);
299: DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);
300: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
301: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
302: switch (dim) {
303: case 2:
304: if (p >= 0) {
305: if (p < nC) {
306: const PetscInt cy = p / nCx;
307: const PetscInt cx = p % nCx;
309: (*cone)[0] = cy * nxF + cx + nC + nV;
310: (*cone)[1] = cx * nyF + cy + nyF + nC + nV + nXF;
311: (*cone)[2] = cy * nxF + cx + nxF + nC + nV;
312: (*cone)[3] = cx * nyF + cy + nC + nV + nXF;
313: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
314: } else if (p < nC + nV) {
315: } else if (p < nC + nV + nXF) {
316: const PetscInt fy = (p - nC + nV) / nxF;
317: const PetscInt fx = (p - nC + nV) % nxF;
319: (*cone)[0] = fy * nVx + fx + nC;
320: (*cone)[1] = fy * nVx + fx + 1 + nC;
321: } else if (p < nC + nV + nXF + nYF) {
322: const PetscInt fx = (p - nC + nV + nXF) / nyF;
323: const PetscInt fy = (p - nC + nV + nXF) % nyF;
325: (*cone)[0] = fy * nVx + fx + nC;
326: (*cone)[1] = fy * nVx + fx + nVx + nC;
327: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " should be in [0, %" PetscInt_FMT ")", p, nC + nV + nXF + nYF + nZF);
328: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %" PetscInt_FMT " is invalid", p);
329: break;
330: case 3:
331: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
332: }
333: return 0;
334: }
336: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
337: {
338: DMGetWorkArray(dm, 6, MPIU_INT, cone);
339: return 0;
340: }
342: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
343: {
344: DM_DA *da = (DM_DA *)dm->data;
345: Vec coordinates;
346: PetscSection section;
347: PetscScalar *coords;
348: PetscReal h[3];
349: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
352: DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
354: h[0] = (xu - xl) / M;
355: h[1] = (yu - yl) / N;
356: h[2] = (zu - zl) / P;
357: DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
358: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
359: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
360: PetscSectionSetNumFields(section, 1);
361: PetscSectionSetFieldComponents(section, 0, dim);
362: PetscSectionSetChart(section, vStart, vEnd);
363: for (v = vStart; v < vEnd; ++v) PetscSectionSetDof(section, v, dim);
364: PetscSectionSetUp(section);
365: PetscSectionGetStorageSize(section, &size);
366: VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
367: PetscObjectSetName((PetscObject)coordinates, "coordinates");
368: VecGetArray(coordinates, &coords);
369: for (k = 0; k < nVz; ++k) {
370: PetscInt ind[3], d, off;
372: ind[0] = 0;
373: ind[1] = 0;
374: ind[2] = k + da->zs;
375: for (j = 0; j < nVy; ++j) {
376: ind[1] = j + da->ys;
377: for (i = 0; i < nVx; ++i) {
378: const PetscInt vertex = (k * nVy + j) * nVx + i + vStart;
380: PetscSectionGetOffset(section, vertex, &off);
381: ind[0] = i + da->xs;
382: for (d = 0; d < dim; ++d) coords[off + d] = h[d] * ind[d];
383: }
384: }
385: }
386: VecRestoreArray(coordinates, &coords);
387: DMSetCoordinateSection(dm, PETSC_DETERMINE, section);
388: DMSetCoordinatesLocal(dm, coordinates);
389: PetscSectionDestroy(§ion);
390: VecDestroy(&coordinates);
391: return 0;
392: }
394: /* ------------------------------------------------------------------- */
396: /*@C
397: DMDAGetArray - Gets a work array for a `DMDA`
399: Input Parameters:
400: + da - information about my local patch
401: - ghosted - do you want arrays for the ghosted or nonghosted patch
403: Output Parameters:
404: . vptr - array data structured
406: Level: advanced
408: Note:
409: The vector values are NOT initialized and may have garbage in them, so you may need
410: to zero them.
412: .seealso: `DM`, `DMDA`, `DMDARestoreArray()`
413: @*/
414: PetscErrorCode DMDAGetArray(DM da, PetscBool ghosted, void *vptr)
415: {
416: PetscInt j, i, xs, ys, xm, ym, zs, zm;
417: char *iarray_start;
418: void **iptr = (void **)vptr;
419: DM_DA *dd = (DM_DA *)da->data;
422: if (ghosted) {
423: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
424: if (dd->arrayghostedin[i]) {
425: *iptr = dd->arrayghostedin[i];
426: iarray_start = (char *)dd->startghostedin[i];
427: dd->arrayghostedin[i] = NULL;
428: dd->startghostedin[i] = NULL;
430: goto done;
431: }
432: }
433: xs = dd->Xs;
434: ys = dd->Ys;
435: zs = dd->Zs;
436: xm = dd->Xe - dd->Xs;
437: ym = dd->Ye - dd->Ys;
438: zm = dd->Ze - dd->Zs;
439: } else {
440: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
441: if (dd->arrayin[i]) {
442: *iptr = dd->arrayin[i];
443: iarray_start = (char *)dd->startin[i];
444: dd->arrayin[i] = NULL;
445: dd->startin[i] = NULL;
447: goto done;
448: }
449: }
450: xs = dd->xs;
451: ys = dd->ys;
452: zs = dd->zs;
453: xm = dd->xe - dd->xs;
454: ym = dd->ye - dd->ys;
455: zm = dd->ze - dd->zs;
456: }
458: switch (da->dim) {
459: case 1: {
460: void *ptr;
462: PetscMalloc(xm * sizeof(PetscScalar), &iarray_start);
464: ptr = (void *)(iarray_start - xs * sizeof(PetscScalar));
465: *iptr = (void *)ptr;
466: break;
467: }
468: case 2: {
469: void **ptr;
471: PetscMalloc((ym + 1) * sizeof(void *) + xm * ym * sizeof(PetscScalar), &iarray_start);
473: ptr = (void **)(iarray_start + xm * ym * sizeof(PetscScalar) - ys * sizeof(void *));
474: for (j = ys; j < ys + ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar) * (xm * (j - ys) - xs);
475: *iptr = (void *)ptr;
476: break;
477: }
478: case 3: {
479: void ***ptr, **bptr;
481: PetscMalloc((zm + 1) * sizeof(void **) + (ym * zm + 1) * sizeof(void *) + xm * ym * zm * sizeof(PetscScalar), &iarray_start);
483: ptr = (void ***)(iarray_start + xm * ym * zm * sizeof(PetscScalar) - zs * sizeof(void *));
484: bptr = (void **)(iarray_start + xm * ym * zm * sizeof(PetscScalar) + zm * sizeof(void **));
485: for (i = zs; i < zs + zm; i++) ptr[i] = bptr + ((i - zs) * ym - ys);
486: for (i = zs; i < zs + zm; i++) {
487: for (j = ys; j < ys + ym; j++) ptr[i][j] = iarray_start + sizeof(PetscScalar) * (xm * ym * (i - zs) + xm * (j - ys) - xs);
488: }
489: *iptr = (void *)ptr;
490: break;
491: }
492: default:
493: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", da->dim);
494: }
496: done:
497: /* add arrays to the checked out list */
498: if (ghosted) {
499: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
500: if (!dd->arrayghostedout[i]) {
501: dd->arrayghostedout[i] = *iptr;
502: dd->startghostedout[i] = iarray_start;
503: break;
504: }
505: }
506: } else {
507: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
508: if (!dd->arrayout[i]) {
509: dd->arrayout[i] = *iptr;
510: dd->startout[i] = iarray_start;
511: break;
512: }
513: }
514: }
515: return 0;
516: }
518: /*@C
519: DMDARestoreArray - Restores an array of derivative types for a `DMDA`
521: Input Parameters:
522: + da - information about my local patch
523: . ghosted - do you want arrays for the ghosted or nonghosted patch
524: - vptr - array data structured
526: Level: advanced
528: .seealso: `DM`, `DMDA`, `DMDAGetArray()`
529: @*/
530: PetscErrorCode DMDARestoreArray(DM da, PetscBool ghosted, void *vptr)
531: {
532: PetscInt i;
533: void **iptr = (void **)vptr, *iarray_start = NULL;
534: DM_DA *dd = (DM_DA *)da->data;
537: if (ghosted) {
538: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
539: if (dd->arrayghostedout[i] == *iptr) {
540: iarray_start = dd->startghostedout[i];
541: dd->arrayghostedout[i] = NULL;
542: dd->startghostedout[i] = NULL;
543: break;
544: }
545: }
546: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
547: if (!dd->arrayghostedin[i]) {
548: dd->arrayghostedin[i] = *iptr;
549: dd->startghostedin[i] = iarray_start;
550: break;
551: }
552: }
553: } else {
554: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
555: if (dd->arrayout[i] == *iptr) {
556: iarray_start = dd->startout[i];
557: dd->arrayout[i] = NULL;
558: dd->startout[i] = NULL;
559: break;
560: }
561: }
562: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
563: if (!dd->arrayin[i]) {
564: dd->arrayin[i] = *iptr;
565: dd->startin[i] = iarray_start;
566: break;
567: }
568: }
569: }
570: return 0;
571: }