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 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: `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: `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 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 - information about my local patch
335: - ghosted - do you want arrays for the ghosted or nonghosted patch
337: Output Parameter:
338: . vptr - array data structured
340: Level: advanced
342: Note:
343: The vector values are NOT initialized and may have garbage in them, so you may need
344: to zero them.
346: .seealso: `DM`, `DMDA`, `DMDARestoreArray()`
347: @*/
348: PetscErrorCode DMDAGetArray(DM da, PetscBool ghosted, void *vptr)
349: {
350: PetscInt j, i, xs, ys, xm, ym, zs, zm;
351: char *iarray_start;
352: void **iptr = (void **)vptr;
353: DM_DA *dd = (DM_DA *)da->data;
355: PetscFunctionBegin;
357: if (ghosted) {
358: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
359: if (dd->arrayghostedin[i]) {
360: *iptr = dd->arrayghostedin[i];
361: iarray_start = (char *)dd->startghostedin[i];
362: dd->arrayghostedin[i] = NULL;
363: dd->startghostedin[i] = NULL;
365: goto done;
366: }
367: }
368: xs = dd->Xs;
369: ys = dd->Ys;
370: zs = dd->Zs;
371: xm = dd->Xe - dd->Xs;
372: ym = dd->Ye - dd->Ys;
373: zm = dd->Ze - dd->Zs;
374: } else {
375: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
376: if (dd->arrayin[i]) {
377: *iptr = dd->arrayin[i];
378: iarray_start = (char *)dd->startin[i];
379: dd->arrayin[i] = NULL;
380: dd->startin[i] = NULL;
382: goto done;
383: }
384: }
385: xs = dd->xs;
386: ys = dd->ys;
387: zs = dd->zs;
388: xm = dd->xe - dd->xs;
389: ym = dd->ye - dd->ys;
390: zm = dd->ze - dd->zs;
391: }
393: switch (da->dim) {
394: case 1: {
395: void *ptr;
397: PetscCall(PetscMalloc(xm * sizeof(PetscScalar), &iarray_start));
399: ptr = (void *)(iarray_start - xs * sizeof(PetscScalar));
400: *iptr = (void *)ptr;
401: break;
402: }
403: case 2: {
404: void **ptr;
406: PetscCall(PetscMalloc((ym + 1) * sizeof(void *) + xm * ym * sizeof(PetscScalar), &iarray_start));
408: ptr = (void **)(iarray_start + xm * ym * sizeof(PetscScalar) - ys * sizeof(void *));
409: for (j = ys; j < ys + ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar) * (xm * (j - ys) - xs);
410: *iptr = (void *)ptr;
411: break;
412: }
413: case 3: {
414: void ***ptr, **bptr;
416: PetscCall(PetscMalloc((zm + 1) * sizeof(void **) + (ym * zm + 1) * sizeof(void *) + xm * ym * zm * sizeof(PetscScalar), &iarray_start));
418: ptr = (void ***)(iarray_start + xm * ym * zm * sizeof(PetscScalar) - zs * sizeof(void *));
419: bptr = (void **)(iarray_start + xm * ym * zm * sizeof(PetscScalar) + zm * sizeof(void **));
420: for (i = zs; i < zs + zm; i++) ptr[i] = bptr + ((i - zs) * ym - ys);
421: for (i = zs; i < zs + zm; i++) {
422: for (j = ys; j < ys + ym; j++) ptr[i][j] = iarray_start + sizeof(PetscScalar) * (xm * ym * (i - zs) + xm * (j - ys) - xs);
423: }
424: *iptr = (void *)ptr;
425: break;
426: }
427: default:
428: SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", da->dim);
429: }
431: done:
432: /* add arrays to the checked out list */
433: if (ghosted) {
434: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
435: if (!dd->arrayghostedout[i]) {
436: dd->arrayghostedout[i] = *iptr;
437: dd->startghostedout[i] = iarray_start;
438: break;
439: }
440: }
441: } else {
442: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
443: if (!dd->arrayout[i]) {
444: dd->arrayout[i] = *iptr;
445: dd->startout[i] = iarray_start;
446: break;
447: }
448: }
449: }
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: /*@C
454: DMDARestoreArray - Restores an array of derivative types for a `DMDA`
456: Input Parameters:
457: + da - information about my local patch
458: . ghosted - do you want arrays for the ghosted or nonghosted patch
459: - vptr - array data structured
461: Level: advanced
463: .seealso: `DM`, `DMDA`, `DMDAGetArray()`
464: @*/
465: PetscErrorCode DMDARestoreArray(DM da, PetscBool ghosted, void *vptr)
466: {
467: PetscInt i;
468: void **iptr = (void **)vptr, *iarray_start = NULL;
469: DM_DA *dd = (DM_DA *)da->data;
471: PetscFunctionBegin;
473: if (ghosted) {
474: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
475: if (dd->arrayghostedout[i] == *iptr) {
476: iarray_start = dd->startghostedout[i];
477: dd->arrayghostedout[i] = NULL;
478: dd->startghostedout[i] = NULL;
479: break;
480: }
481: }
482: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
483: if (!dd->arrayghostedin[i]) {
484: dd->arrayghostedin[i] = *iptr;
485: dd->startghostedin[i] = iarray_start;
486: break;
487: }
488: }
489: } else {
490: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
491: if (dd->arrayout[i] == *iptr) {
492: iarray_start = dd->startout[i];
493: dd->arrayout[i] = NULL;
494: dd->startout[i] = NULL;
495: break;
496: }
497: }
498: for (i = 0; i < DMDA_MAX_WORK_ARRAYS; i++) {
499: if (!dd->arrayin[i]) {
500: dd->arrayin[i] = *iptr;
501: dd->startin[i] = iarray_start;
502: break;
503: }
504: }
505: }
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }