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_ENGINE)
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_ENGINE)
61: if (dd->w == 1 && da->dim == 2) {
62: PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);
63: }
64: #endif
65: return 0;
66: }
68: /*@
69: DMDAGetNumCells - Get the number of cells in the local piece of the DMDA. This includes ghost cells.
71: Input Parameter:
72: . dm - The DM object
74: Output Parameters:
75: + numCellsX - The number of local cells in the x-direction
76: . numCellsY - The number of local cells in the y-direction
77: . numCellsZ - The number of local cells in the z-direction
78: - numCells - The number of local cells
80: Level: developer
82: .seealso: DMDAGetCellPoint()
83: @*/
84: PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells)
85: {
86: DM_DA *da = (DM_DA*) dm->data;
87: const PetscInt dim = dm->dim;
88: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
89: const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
92: if (numCellsX) {
94: *numCellsX = mx;
95: }
96: if (numCellsY) {
98: *numCellsY = my;
99: }
100: if (numCellsZ) {
102: *numCellsZ = mz;
103: }
104: if (numCells) {
106: *numCells = nC;
107: }
108: return 0;
109: }
111: /*@
112: DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the DMDA
114: Input Parameters:
115: + dm - The DM object
116: - i,j,k - The global indices for the cell
118: Output Parameters:
119: . point - The local DM point
121: Level: developer
123: .seealso: DMDAGetNumCells()
124: @*/
125: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
126: {
127: const PetscInt dim = dm->dim;
128: DMDALocalInfo info;
132: DMDAGetLocalInfo(dm, &info);
136: *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0);
137: return 0;
138: }
140: PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices)
141: {
142: DM_DA *da = (DM_DA*) dm->data;
143: const PetscInt dim = dm->dim;
144: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
145: const PetscInt nVx = mx+1;
146: const PetscInt nVy = dim > 1 ? (my+1) : 1;
147: const PetscInt nVz = dim > 2 ? (mz+1) : 1;
148: const PetscInt nV = nVx*nVy*nVz;
150: if (numVerticesX) {
152: *numVerticesX = nVx;
153: }
154: if (numVerticesY) {
156: *numVerticesY = nVy;
157: }
158: if (numVerticesZ) {
160: *numVerticesZ = nVz;
161: }
162: if (numVertices) {
164: *numVertices = nV;
165: }
166: return 0;
167: }
169: PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces)
170: {
171: DM_DA *da = (DM_DA*) dm->data;
172: const PetscInt dim = dm->dim;
173: const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs;
174: const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1);
175: const PetscInt nXF = (mx+1)*nxF;
176: const PetscInt nyF = mx*(dim > 2 ? mz : 1);
177: const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0;
178: const PetscInt nzF = mx*(dim > 1 ? my : 0);
179: const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0;
181: if (numXFacesX) {
183: *numXFacesX = nxF;
184: }
185: if (numXFaces) {
187: *numXFaces = nXF;
188: }
189: if (numYFacesY) {
191: *numYFacesY = nyF;
192: }
193: if (numYFaces) {
195: *numYFaces = nYF;
196: }
197: if (numZFacesZ) {
199: *numZFacesZ = nzF;
200: }
201: if (numZFaces) {
203: *numZFaces = nZF;
204: }
205: return 0;
206: }
208: PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd)
209: {
210: const PetscInt dim = dm->dim;
211: PetscInt nC, nV, nXF, nYF, nZF;
215: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
216: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
217: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
218: if (height == 0) {
219: /* Cells */
220: if (pStart) *pStart = 0;
221: if (pEnd) *pEnd = nC;
222: } else if (height == 1) {
223: /* Faces */
224: if (pStart) *pStart = nC+nV;
225: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
226: } else if (height == dim) {
227: /* Vertices */
228: if (pStart) *pStart = nC;
229: if (pEnd) *pEnd = nC+nV;
230: } else if (height < 0) {
231: /* All points */
232: if (pStart) *pStart = 0;
233: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
234: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
235: return 0;
236: }
238: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
239: {
240: const PetscInt dim = dm->dim;
241: PetscInt nC, nV, nXF, nYF, nZF;
245: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
246: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
247: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
248: if (depth == dim) {
249: /* Cells */
250: if (pStart) *pStart = 0;
251: if (pEnd) *pEnd = nC;
252: } else if (depth == dim-1) {
253: /* Faces */
254: if (pStart) *pStart = nC+nV;
255: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
256: } else if (depth == 0) {
257: /* Vertices */
258: if (pStart) *pStart = nC;
259: if (pEnd) *pEnd = nC+nV;
260: } else if (depth < 0) {
261: /* All points */
262: if (pStart) *pStart = 0;
263: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
264: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
265: return 0;
266: }
268: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
269: {
270: const PetscInt dim = dm->dim;
271: PetscInt nC, nV, nXF, nYF, nZF;
273: *coneSize = 0;
274: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
275: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
276: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
277: switch (dim) {
278: case 2:
279: if (p >= 0) {
280: if (p < nC) {
281: *coneSize = 4;
282: } else if (p < nC+nV) {
283: *coneSize = 0;
284: } else if (p < nC+nV+nXF+nYF+nZF) {
285: *coneSize = 2;
286: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
287: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
288: break;
289: case 3:
290: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
291: }
292: return 0;
293: }
295: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
296: {
297: const PetscInt dim = dm->dim;
298: PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
300: if (!cone) DMGetWorkArray(dm, 6, MPIU_INT, cone);
301: DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);
302: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
303: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
304: switch (dim) {
305: case 2:
306: if (p >= 0) {
307: if (p < nC) {
308: const PetscInt cy = p / nCx;
309: const PetscInt cx = p % nCx;
311: (*cone)[0] = cy*nxF + cx + nC+nV;
312: (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
313: (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
314: (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
315: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
316: } else if (p < nC+nV) {
317: } else if (p < nC+nV+nXF) {
318: const PetscInt fy = (p - nC+nV) / nxF;
319: const PetscInt fx = (p - nC+nV) % nxF;
321: (*cone)[0] = fy*nVx + fx + nC;
322: (*cone)[1] = fy*nVx + fx + 1 + nC;
323: } else if (p < nC+nV+nXF+nYF) {
324: const PetscInt fx = (p - nC+nV+nXF) / nyF;
325: const PetscInt fy = (p - nC+nV+nXF) % nyF;
327: (*cone)[0] = fy*nVx + fx + nC;
328: (*cone)[1] = fy*nVx + fx + nVx + nC;
329: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
330: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
331: break;
332: case 3:
333: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
334: }
335: return 0;
336: }
338: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
339: {
340: DMGetWorkArray(dm, 6, MPIU_INT, cone);
341: return 0;
342: }
344: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
345: {
346: DM_DA *da = (DM_DA *) dm->data;
347: Vec coordinates;
348: PetscSection section;
349: PetscScalar *coords;
350: PetscReal h[3];
351: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
354: DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
356: h[0] = (xu - xl)/M;
357: h[1] = (yu - yl)/N;
358: h[2] = (zu - zl)/P;
359: DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
360: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
361: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);
362: PetscSectionSetNumFields(section, 1);
363: PetscSectionSetFieldComponents(section, 0, dim);
364: PetscSectionSetChart(section, vStart, vEnd);
365: for (v = vStart; v < vEnd; ++v) {
366: PetscSectionSetDof(section, v, dim);
367: }
368: PetscSectionSetUp(section);
369: PetscSectionGetStorageSize(section, &size);
370: VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
371: PetscObjectSetName((PetscObject)coordinates,"coordinates");
372: VecGetArray(coordinates, &coords);
373: for (k = 0; k < nVz; ++k) {
374: PetscInt ind[3], d, off;
376: ind[0] = 0;
377: ind[1] = 0;
378: ind[2] = k + da->zs;
379: for (j = 0; j < nVy; ++j) {
380: ind[1] = j + da->ys;
381: for (i = 0; i < nVx; ++i) {
382: const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;
384: PetscSectionGetOffset(section, vertex, &off);
385: ind[0] = i + da->xs;
386: for (d = 0; d < dim; ++d) {
387: coords[off+d] = h[d]*ind[d];
388: }
389: }
390: }
391: }
392: VecRestoreArray(coordinates, &coords);
393: DMSetCoordinateSection(dm, PETSC_DETERMINE, section);
394: DMSetCoordinatesLocal(dm, coordinates);
395: PetscSectionDestroy(§ion);
396: VecDestroy(&coordinates);
397: return 0;
398: }
400: /* ------------------------------------------------------------------- */
402: /*@C
403: DMDAGetArray - Gets a work array for a DMDA
405: Input Parameters:
406: + da - information about my local patch
407: - ghosted - do you want arrays for the ghosted or nonghosted patch
409: Output Parameters:
410: . vptr - array data structured
412: Note: The vector values are NOT initialized and may have garbage in them, so you may need
413: to zero them.
415: Level: advanced
417: .seealso: DMDARestoreArray()
419: @*/
420: PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
421: {
422: PetscInt j,i,xs,ys,xm,ym,zs,zm;
423: char *iarray_start;
424: void **iptr = (void**)vptr;
425: DM_DA *dd = (DM_DA*)da->data;
428: if (ghosted) {
429: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
430: if (dd->arrayghostedin[i]) {
431: *iptr = dd->arrayghostedin[i];
432: iarray_start = (char*)dd->startghostedin[i];
433: dd->arrayghostedin[i] = NULL;
434: dd->startghostedin[i] = NULL;
436: goto done;
437: }
438: }
439: xs = dd->Xs;
440: ys = dd->Ys;
441: zs = dd->Zs;
442: xm = dd->Xe-dd->Xs;
443: ym = dd->Ye-dd->Ys;
444: zm = dd->Ze-dd->Zs;
445: } else {
446: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
447: if (dd->arrayin[i]) {
448: *iptr = dd->arrayin[i];
449: iarray_start = (char*)dd->startin[i];
450: dd->arrayin[i] = NULL;
451: dd->startin[i] = NULL;
453: goto done;
454: }
455: }
456: xs = dd->xs;
457: ys = dd->ys;
458: zs = dd->zs;
459: xm = dd->xe-dd->xs;
460: ym = dd->ye-dd->ys;
461: zm = dd->ze-dd->zs;
462: }
464: switch (da->dim) {
465: case 1: {
466: void *ptr;
468: PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);
470: ptr = (void*)(iarray_start - xs*sizeof(PetscScalar));
471: *iptr = (void*)ptr;
472: break;
473: }
474: case 2: {
475: void **ptr;
477: PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);
479: ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
480: for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
481: *iptr = (void*)ptr;
482: break;
483: }
484: case 3: {
485: void ***ptr,**bptr;
487: PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);
489: ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
490: bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
491: for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
492: for (i=zs; i<zs+zm; i++) {
493: for (j=ys; j<ys+ym; j++) {
494: ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
495: }
496: }
497: *iptr = (void*)ptr;
498: break;
499: }
500: default:
501: SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
502: }
504: done:
505: /* add arrays to the checked out list */
506: if (ghosted) {
507: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
508: if (!dd->arrayghostedout[i]) {
509: dd->arrayghostedout[i] = *iptr;
510: dd->startghostedout[i] = iarray_start;
511: break;
512: }
513: }
514: } else {
515: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
516: if (!dd->arrayout[i]) {
517: dd->arrayout[i] = *iptr;
518: dd->startout[i] = iarray_start;
519: break;
520: }
521: }
522: }
523: return 0;
524: }
526: /*@C
527: DMDARestoreArray - Restores an array of derivative types for a DMDA
529: Input Parameters:
530: + da - information about my local patch
531: . ghosted - do you want arrays for the ghosted or nonghosted patch
532: - vptr - array data structured
534: Level: advanced
536: .seealso: DMDAGetArray()
538: @*/
539: PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
540: {
541: PetscInt i;
542: void **iptr = (void**)vptr,*iarray_start = NULL;
543: DM_DA *dd = (DM_DA*)da->data;
546: if (ghosted) {
547: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
548: if (dd->arrayghostedout[i] == *iptr) {
549: iarray_start = dd->startghostedout[i];
550: dd->arrayghostedout[i] = NULL;
551: dd->startghostedout[i] = NULL;
552: break;
553: }
554: }
555: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
556: if (!dd->arrayghostedin[i]) {
557: dd->arrayghostedin[i] = *iptr;
558: dd->startghostedin[i] = iarray_start;
559: break;
560: }
561: }
562: } else {
563: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
564: if (dd->arrayout[i] == *iptr) {
565: iarray_start = dd->startout[i];
566: dd->arrayout[i] = NULL;
567: dd->startout[i] = NULL;
568: break;
569: }
570: }
571: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
572: if (!dd->arrayin[i]) {
573: dd->arrayin[i] = *iptr;
574: dd->startin[i] = iarray_start;
575: break;
576: }
577: }
578: }
579: return 0;
580: }