Actual source code: dalocal.c
1: /*
2: Code for manipulating distributed regular arrays in parallel.
3: */
5: #include <petsc/private/dmdaimpl.h>
6: #include <petscbt.h>
7: #include <petscsf.h>
8: #include <petscds.h>
9: #include <petscfe.h>
11: /*
12: This allows the DMDA vectors to properly tell MATLAB their dimensions
13: */
14: #if defined(PETSC_HAVE_MATLAB)
15: #include <engine.h> /* MATLAB include file */
16: #include <mex.h> /* MATLAB include file */
17: static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj, void *mengine)
18: {
19: PetscInt n, m;
20: Vec vec = (Vec)obj;
21: PetscScalar *array;
22: mxArray *mat;
23: DM da;
25: PetscFunctionBegin;
26: PetscCall(VecGetDM(vec, &da));
27: PetscCheck(da, PetscObjectComm((PetscObject)vec), PETSC_ERR_ARG_WRONGSTATE, "Vector not associated with a DMDA");
28: PetscCall(DMDAGetGhostCorners(da, 0, 0, 0, &m, &n, 0));
30: PetscCall(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: PetscCall(PetscArraycpy(mxGetPr(mat), array, n * m));
37: PetscCall(PetscObjectName(obj));
38: engPutVariable((Engine *)mengine, obj->name, mat);
40: PetscCall(VecRestoreArray(vec, &array));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
43: #endif
45: PetscErrorCode DMCreateLocalVector_DA(DM da, Vec *g)
46: {
47: DM_DA *dd = (DM_DA *)da->data;
49: PetscFunctionBegin;
51: PetscAssertPointer(g, 2);
52: PetscCall(VecCreate(PETSC_COMM_SELF, g));
53: PetscCall(VecSetSizes(*g, dd->nlocal, PETSC_DETERMINE));
54: PetscCall(VecSetBlockSize(*g, dd->w));
55: PetscCall(VecSetType(*g, da->vectype));
56: if (dd->nlocal < da->bind_below) {
57: PetscCall(VecSetBindingPropagates(*g, PETSC_TRUE));
58: PetscCall(VecBindToCPU(*g, PETSC_TRUE));
59: }
60: PetscCall(VecSetDM(*g, da));
61: #if defined(PETSC_HAVE_MATLAB)
62: if (dd->w == 1 && da->dim == 2) PetscCall(PetscObjectComposeFunction((PetscObject)*g, "PetscMatlabEnginePut_C", VecMatlabEnginePut_DA2d));
63: #endif
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*@
68: DMDAGetNumCells - Get the number of cells (or vertices) in the local piece of the `DMDA`. This includes ghost cells.
70: Input Parameter:
71: . dm - The `DMDA` object
73: Output Parameters:
74: + numCellsX - The number of local cells in the x-direction
75: . numCellsY - The number of local cells in the y-direction
76: . numCellsZ - The number of local cells in the z-direction
77: - numCells - The number of local cells
79: Level: developer
81: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetCellPoint()`
82: @*/
83: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
84: {
85: DM_DA *da = (DM_DA *)dm->data;
86: const PetscInt dim = dm->dim;
87: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
88: const PetscInt nC = (mx) * (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
90: PetscFunctionBegin;
92: if (numCellsX) {
93: PetscAssertPointer(numCellsX, 2);
94: *numCellsX = mx;
95: }
96: if (numCellsY) {
97: PetscAssertPointer(numCellsY, 3);
98: *numCellsY = my;
99: }
100: if (numCellsZ) {
101: PetscAssertPointer(numCellsZ, 4);
102: *numCellsZ = mz;
103: }
104: if (numCells) {
105: PetscAssertPointer(numCells, 5);
106: *numCells = nC;
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: /*@
112: DMDAGetCellPoint - Get the `DM` point corresponding to the tuple (i, j, k) in the `DMDA`
114: Input Parameters:
115: + dm - The `DMDA` object
116: . i - The global x index for the cell
117: . j - The global y index for the cell
118: - k - The global z index for the cell
120: Output Parameter:
121: . point - The local `DM` point
123: Level: developer
125: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetNumCells()`
126: @*/
127: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
128: {
129: const PetscInt dim = dm->dim;
130: DMDALocalInfo info;
132: PetscFunctionBegin;
134: PetscAssertPointer(point, 5);
135: PetscCall(DMDAGetLocalInfo(dm, &info));
136: if (dim > 0) PetscCheck(!(i < info.gxs) && !(i >= info.gxs + info.gxm), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "X index %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", i, info.gxs, info.gxs + info.gxm);
137: if (dim > 1) PetscCheck(!(j < info.gys) && !(j >= info.gys + info.gym), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", j, info.gys, info.gys + info.gym);
138: if (dim > 2) PetscCheck(!(k < info.gzs) && !(k >= info.gzs + info.gzm), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %" PetscInt_FMT " not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", k, info.gzs, info.gzs + info.gzm);
139: *point = i + (dim > 1 ? (j + (dim > 2 ? k * info.gym : 0)) * info.gxm : 0);
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
144: {
145: DM_DA *da = (DM_DA *)dm->data;
146: const PetscInt dim = dm->dim;
147: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
148: const PetscInt nVx = mx + 1;
149: const PetscInt nVy = dim > 1 ? (my + 1) : 1;
150: const PetscInt nVz = dim > 2 ? (mz + 1) : 1;
151: const PetscInt nV = nVx * nVy * nVz;
153: PetscFunctionBegin;
154: if (numVerticesX) {
155: PetscAssertPointer(numVerticesX, 2);
156: *numVerticesX = nVx;
157: }
158: if (numVerticesY) {
159: PetscAssertPointer(numVerticesY, 3);
160: *numVerticesY = nVy;
161: }
162: if (numVerticesZ) {
163: PetscAssertPointer(numVerticesZ, 4);
164: *numVerticesZ = nVz;
165: }
166: if (numVertices) {
167: PetscAssertPointer(numVertices, 5);
168: *numVertices = nV;
169: }
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
174: {
175: DM_DA *da = (DM_DA *)dm->data;
176: const PetscInt dim = dm->dim;
177: const PetscInt mx = (da->Xe - da->Xs) / da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
178: const PetscInt nxF = (dim > 1 ? (my) * (dim > 2 ? (mz) : 1) : 1);
179: const PetscInt nXF = (mx + 1) * nxF;
180: const PetscInt nyF = mx * (dim > 2 ? mz : 1);
181: const PetscInt nYF = dim > 1 ? (my + 1) * nyF : 0;
182: const PetscInt nzF = mx * (dim > 1 ? my : 0);
183: const PetscInt nZF = dim > 2 ? (mz + 1) * nzF : 0;
185: PetscFunctionBegin;
186: if (numXFacesX) {
187: PetscAssertPointer(numXFacesX, 2);
188: *numXFacesX = nxF;
189: }
190: if (numXFaces) {
191: PetscAssertPointer(numXFaces, 3);
192: *numXFaces = nXF;
193: }
194: if (numYFacesY) {
195: PetscAssertPointer(numYFacesY, 4);
196: *numYFacesY = nyF;
197: }
198: if (numYFaces) {
199: PetscAssertPointer(numYFaces, 5);
200: *numYFaces = nYF;
201: }
202: if (numZFacesZ) {
203: PetscAssertPointer(numZFacesZ, 6);
204: *numZFacesZ = nzF;
205: }
206: if (numZFaces) {
207: PetscAssertPointer(numZFaces, 7);
208: *numZFaces = nZF;
209: }
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
214: {
215: const PetscInt dim = dm->dim;
216: PetscInt nC, nV, nXF, nYF, nZF;
218: PetscFunctionBegin;
219: if (pStart) PetscAssertPointer(pStart, 3);
220: if (pEnd) PetscAssertPointer(pEnd, 4);
221: PetscCall(DMDAGetNumCells(dm, NULL, NULL, NULL, &nC));
222: PetscCall(DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV));
223: PetscCall(DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF));
224: if (height == 0) {
225: /* Cells */
226: if (pStart) *pStart = 0;
227: if (pEnd) *pEnd = nC;
228: } else if (height == 1) {
229: /* Faces */
230: if (pStart) *pStart = nC + nV;
231: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
232: } else if (height == dim) {
233: /* Vertices */
234: if (pStart) *pStart = nC;
235: if (pEnd) *pEnd = nC + nV;
236: } else if (height < 0) {
237: /* All points */
238: if (pStart) *pStart = 0;
239: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
240: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %" PetscInt_FMT " in the DA", height);
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
245: {
246: const PetscInt dim = dm->dim;
247: PetscInt nC, nV, nXF, nYF, nZF;
249: PetscFunctionBegin;
250: if (pStart) PetscAssertPointer(pStart, 3);
251: if (pEnd) PetscAssertPointer(pEnd, 4);
252: PetscCall(DMDAGetNumCells(dm, NULL, NULL, NULL, &nC));
253: PetscCall(DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV));
254: PetscCall(DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF));
255: if (depth == dim) {
256: /* Cells */
257: if (pStart) *pStart = 0;
258: if (pEnd) *pEnd = nC;
259: } else if (depth == dim - 1) {
260: /* Faces */
261: if (pStart) *pStart = nC + nV;
262: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
263: } else if (depth == 0) {
264: /* Vertices */
265: if (pStart) *pStart = nC;
266: if (pEnd) *pEnd = nC + nV;
267: } else if (depth < 0) {
268: /* All points */
269: if (pStart) *pStart = 0;
270: if (pEnd) *pEnd = nC + nV + nXF + nYF + nZF;
271: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %" PetscInt_FMT " in the DA", depth);
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
276: {
277: DM_DA *da = (DM_DA *)dm->data;
278: Vec coordinates;
279: PetscSection section;
280: PetscScalar *coords;
281: PetscReal h[3];
282: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
284: PetscFunctionBegin;
286: PetscCall(DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
287: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The following code only works for dim <= 3");
288: h[0] = (xu - xl) / M;
289: h[1] = (yu - yl) / N;
290: h[2] = (zu - zl) / P;
291: PetscCall(DMDAGetDepthStratum(dm, 0, &vStart, &vEnd));
292: PetscCall(DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV));
293: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
294: PetscCall(PetscSectionSetNumFields(section, 1));
295: PetscCall(PetscSectionSetFieldComponents(section, 0, dim));
296: PetscCall(PetscSectionSetChart(section, vStart, vEnd));
297: for (v = vStart; v < vEnd; ++v) PetscCall(PetscSectionSetDof(section, v, dim));
298: PetscCall(PetscSectionSetUp(section));
299: PetscCall(PetscSectionGetStorageSize(section, &size));
300: PetscCall(VecCreateSeq(PETSC_COMM_SELF, size, &coordinates));
301: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
302: PetscCall(VecGetArray(coordinates, &coords));
303: for (k = 0; k < nVz; ++k) {
304: PetscInt ind[3], d, off;
306: ind[0] = 0;
307: ind[1] = 0;
308: ind[2] = k + da->zs;
309: for (j = 0; j < nVy; ++j) {
310: ind[1] = j + da->ys;
311: for (i = 0; i < nVx; ++i) {
312: const PetscInt vertex = (k * nVy + j) * nVx + i + vStart;
314: PetscCall(PetscSectionGetOffset(section, vertex, &off));
315: ind[0] = i + da->xs;
316: for (d = 0; d < dim; ++d) coords[off + d] = h[d] * ind[d];
317: }
318: }
319: }
320: PetscCall(VecRestoreArray(coordinates, &coords));
321: PetscCall(DMSetCoordinateSection(dm, PETSC_DETERMINE, section));
322: PetscCall(DMSetCoordinatesLocal(dm, coordinates));
323: PetscCall(PetscSectionDestroy(§ion));
324: PetscCall(VecDestroy(&coordinates));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /* ------------------------------------------------------------------- */
330: /*@C
331: DMDAGetArray - Gets a work array for a `DMDA`
333: Input Parameters:
334: + da - a `DMDA`
335: - ghosted - do you want arrays for the ghosted or nonghosted patch
337: Output Parameter:
338: . vptr - array data structured
340: Level: advanced
342: Notes:
343: The vector values are NOT initialized and may have garbage in them, so you may need
344: to zero them.
346: Use `DMDARestoreArray()` to return the array
348: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDARestoreArray()`
349: @*/
350: PetscErrorCode DMDAGetArray(DM da, PetscBool ghosted, void *vptr)
351: {
352: PetscInt j, i, xs, ys, xm, ym, zs, zm;
353: char *iarray_start;
354: void **iptr = (void **)vptr;
355: DM_DA *dd = (DM_DA *)da->data;
357: PetscFunctionBegin;
359: if (ghosted) {
360: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
361: if (dd->arrayghostedin[i]) {
362: *iptr = dd->arrayghostedin[i];
363: iarray_start = (char *)dd->startghostedin[i];
364: dd->arrayghostedin[i] = NULL;
365: dd->startghostedin[i] = NULL;
367: goto done;
368: }
369: }
370: xs = dd->Xs;
371: ys = dd->Ys;
372: zs = dd->Zs;
373: xm = dd->Xe - dd->Xs;
374: ym = dd->Ye - dd->Ys;
375: zm = dd->Ze - dd->Zs;
376: } else {
377: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
378: if (dd->arrayin[i]) {
379: *iptr = dd->arrayin[i];
380: iarray_start = (char *)dd->startin[i];
381: dd->arrayin[i] = NULL;
382: dd->startin[i] = NULL;
384: goto done;
385: }
386: }
387: xs = dd->xs;
388: ys = dd->ys;
389: zs = dd->zs;
390: xm = dd->xe - dd->xs;
391: ym = dd->ye - dd->ys;
392: zm = dd->ze - dd->zs;
393: }
395: switch (da->dim) {
396: case 1: {
397: void *ptr;
399: PetscCall(PetscMalloc(xm * sizeof(PetscScalar), &iarray_start));
401: ptr = (void *)(iarray_start - xs * sizeof(PetscScalar));
402: *iptr = (void *)ptr;
403: break;
404: }
405: case 2: {
406: void **ptr;
408: PetscCall(PetscMalloc((ym + 1) * sizeof(void *) + xm * ym * sizeof(PetscScalar), &iarray_start));
410: ptr = (void **)(iarray_start + xm * ym * sizeof(PetscScalar) - ys * sizeof(void *));
411: for (j = ys; j < ys + ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar) * (xm * (j - ys) - xs);
412: *iptr = (void *)ptr;
413: break;
414: }
415: case 3: {
416: void ***ptr, **bptr;
418: PetscCall(PetscMalloc((zm + 1) * sizeof(void **) + (ym * zm + 1) * sizeof(void *) + xm * ym * zm * sizeof(PetscScalar), &iarray_start));
420: ptr = (void ***)(iarray_start + xm * ym * zm * sizeof(PetscScalar) - zs * sizeof(void *));
421: bptr = (void **)(iarray_start + xm * ym * zm * sizeof(PetscScalar) + zm * sizeof(void **));
422: for (i = zs; i < zs + zm; i++) ptr[i] = bptr + ((i - zs) * ym - ys);
423: for (i = zs; i < zs + zm; i++) {
424: for (j = ys; j < ys + ym; j++) ptr[i][j] = iarray_start + sizeof(PetscScalar) * (xm * ym * (i - zs) + xm * (j - ys) - xs);
425: }
426: *iptr = (void *)ptr;
427: break;
428: }
429: default:
430: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", da->dim);
431: }
433: done:
434: /* add arrays to the checked out list */
435: if (ghosted) {
436: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
437: if (!dd->arrayghostedout[i]) {
438: dd->arrayghostedout[i] = *iptr;
439: dd->startghostedout[i] = iarray_start;
440: break;
441: }
442: }
443: } else {
444: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
445: if (!dd->arrayout[i]) {
446: dd->arrayout[i] = *iptr;
447: dd->startout[i] = iarray_start;
448: break;
449: }
450: }
451: }
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: /*@C
456: DMDARestoreArray - Restores an array for a `DMDA` obtained with `DMDAGetArray()`
458: Input Parameters:
459: + da - information about my local patch
460: . ghosted - do you want arrays for the ghosted or nonghosted patch
461: - vptr - array data structured
463: Level: advanced
465: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDAGetArray()`
466: @*/
467: PetscErrorCode DMDARestoreArray(DM da, PetscBool ghosted, void *vptr)
468: {
469: PetscInt i;
470: void **iptr = (void **)vptr, *iarray_start = NULL;
471: DM_DA *dd = (DM_DA *)da->data;
473: PetscFunctionBegin;
475: if (ghosted) {
476: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
477: if (dd->arrayghostedout[i] == *iptr) {
478: iarray_start = dd->startghostedout[i];
479: dd->arrayghostedout[i] = NULL;
480: dd->startghostedout[i] = NULL;
481: break;
482: }
483: }
484: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
485: if (!dd->arrayghostedin[i]) {
486: dd->arrayghostedin[i] = *iptr;
487: dd->startghostedin[i] = iarray_start;
488: break;
489: }
490: }
491: } else {
492: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
493: if (dd->arrayout[i] == *iptr) {
494: iarray_start = dd->startout[i];
495: dd->arrayout[i] = NULL;
496: dd->startout[i] = NULL;
497: break;
498: }
499: }
500: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
501: if (!dd->arrayin[i]) {
502: dd->arrayin[i] = *iptr;
503: dd->startin[i] = iarray_start;
504: break;
505: }
506: }
507: }
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }