Actual source code: dalocal.c
petsc-3.13.6 2020-09-29
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: {
21: PetscInt n,m;
22: Vec vec = (Vec)obj;
23: PetscScalar *array;
24: mxArray *mat;
25: DM da;
28: VecGetDM(vec, &da);
29: if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA");
30: DMDAGetGhostCorners(da,0,0,0,&m,&n,0);
32: VecGetArray(vec,&array);
33: #if !defined(PETSC_USE_COMPLEX)
34: mat = mxCreateDoubleMatrix(m,n,mxREAL);
35: #else
36: mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX);
37: #endif
38: PetscArraycpy(mxGetPr(mat),array,n*m);
39: PetscObjectName(obj);
40: engPutVariable((Engine*)mengine,obj->name,mat);
42: VecRestoreArray(vec,&array);
43: return(0);
44: }
45: #endif
47: PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g)
48: {
50: DM_DA *dd = (DM_DA*)da->data;
55: VecCreate(PETSC_COMM_SELF,g);
56: VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);
57: VecSetBlockSize(*g,dd->w);
58: VecSetType(*g,da->vectype);
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);
93: if (numCellsX) {
95: *numCellsX = mx;
96: }
97: if (numCellsY) {
99: *numCellsY = my;
100: }
101: if (numCellsZ) {
103: *numCellsZ = mz;
104: }
105: if (numCells) {
107: *numCells = nC;
108: }
109: return(0);
110: }
112: /*@
113: DMDAGetCellPoint - Get the DM point corresponding to the tuple (i, j, k) in the DMDA
115: Input Parameters:
116: + dm - The DM object
117: - i,j,k - The global indices for the cell
119: Output Parameters:
120: . point - The local DM point
122: Level: developer
124: .seealso: DMDAGetNumCells()
125: @*/
126: PetscErrorCode DMDAGetCellPoint(DM dm, PetscInt i, PetscInt j, PetscInt k, PetscInt *point)
127: {
128: const PetscInt dim = dm->dim;
129: DMDALocalInfo info;
135: DMDAGetLocalInfo(dm, &info);
136: if (dim > 0) {if ((i < info.gxs) || (i >= info.gxs+info.gxm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "X index %d not in [%d, %d)", i, info.gxs, info.gxs+info.gxm);}
137: if (dim > 1) {if ((j < info.gys) || (j >= info.gys+info.gym)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Y index %d not in [%d, %d)", j, info.gys, info.gys+info.gym);}
138: if (dim > 2) {if ((k < info.gzs) || (k >= info.gzs+info.gzm)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Z index %d not in [%d, %d)", k, info.gzs, info.gzs+info.gzm);}
139: *point = i + (dim > 1 ? (j + (dim > 2 ? k*info.gym : 0))*info.gxm : 0);
140: return(0);
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;
154: if (numVerticesX) {
156: *numVerticesX = nVx;
157: }
158: if (numVerticesY) {
160: *numVerticesY = nVy;
161: }
162: if (numVerticesZ) {
164: *numVerticesZ = nVz;
165: }
166: if (numVertices) {
168: *numVertices = nV;
169: }
170: return(0);
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;
186: if (numXFacesX) {
188: *numXFacesX = nxF;
189: }
190: if (numXFaces) {
192: *numXFaces = nXF;
193: }
194: if (numYFacesY) {
196: *numYFacesY = nyF;
197: }
198: if (numYFaces) {
200: *numYFaces = nYF;
201: }
202: if (numZFacesZ) {
204: *numZFacesZ = nzF;
205: }
206: if (numZFaces) {
208: *numZFaces = nZF;
209: }
210: return(0);
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;
222: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
223: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
224: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
225: if (height == 0) {
226: /* Cells */
227: if (pStart) *pStart = 0;
228: if (pEnd) *pEnd = nC;
229: } else if (height == 1) {
230: /* Faces */
231: if (pStart) *pStart = nC+nV;
232: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
233: } else if (height == dim) {
234: /* Vertices */
235: if (pStart) *pStart = nC;
236: if (pEnd) *pEnd = nC+nV;
237: } else if (height < 0) {
238: /* All points */
239: if (pStart) *pStart = 0;
240: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
241: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height);
242: return(0);
243: }
245: PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd)
246: {
247: const PetscInt dim = dm->dim;
248: PetscInt nC, nV, nXF, nYF, nZF;
254: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
255: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
256: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
257: if (depth == dim) {
258: /* Cells */
259: if (pStart) *pStart = 0;
260: if (pEnd) *pEnd = nC;
261: } else if (depth == dim-1) {
262: /* Faces */
263: if (pStart) *pStart = nC+nV;
264: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
265: } else if (depth == 0) {
266: /* Vertices */
267: if (pStart) *pStart = nC;
268: if (pEnd) *pEnd = nC+nV;
269: } else if (depth < 0) {
270: /* All points */
271: if (pStart) *pStart = 0;
272: if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF;
273: } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth);
274: return(0);
275: }
277: PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize)
278: {
279: const PetscInt dim = dm->dim;
280: PetscInt nC, nV, nXF, nYF, nZF;
284: *coneSize = 0;
285: DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);
286: DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);
287: DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);
288: switch (dim) {
289: case 2:
290: if (p >= 0) {
291: if (p < nC) {
292: *coneSize = 4;
293: } else if (p < nC+nV) {
294: *coneSize = 0;
295: } else if (p < nC+nV+nXF+nYF+nZF) {
296: *coneSize = 2;
297: } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
298: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
299: break;
300: case 3:
301: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
302: break;
303: }
304: return(0);
305: }
307: PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[])
308: {
309: const PetscInt dim = dm->dim;
310: PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF;
314: if (!cone) {DMGetWorkArray(dm, 6, MPIU_INT, cone);}
315: DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);
316: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
317: DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);
318: switch (dim) {
319: case 2:
320: if (p >= 0) {
321: if (p < nC) {
322: const PetscInt cy = p / nCx;
323: const PetscInt cx = p % nCx;
325: (*cone)[0] = cy*nxF + cx + nC+nV;
326: (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF;
327: (*cone)[2] = cy*nxF + cx + nxF + nC+nV;
328: (*cone)[3] = cx*nyF + cy + nC+nV+nXF;
329: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones");
330: } else if (p < nC+nV) {
331: } else if (p < nC+nV+nXF) {
332: const PetscInt fy = (p - nC+nV) / nxF;
333: const PetscInt fx = (p - nC+nV) % nxF;
335: (*cone)[0] = fy*nVx + fx + nC;
336: (*cone)[1] = fy*nVx + fx + 1 + nC;
337: } else if (p < nC+nV+nXF+nYF) {
338: const PetscInt fx = (p - nC+nV+nXF) / nyF;
339: const PetscInt fy = (p - nC+nV+nXF) % nyF;
341: (*cone)[0] = fy*nVx + fx + nC;
342: (*cone)[1] = fy*nVx + fx + nVx + nC;
343: } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF);
344: } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p);
345: break;
346: case 3:
347: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D");
348: break;
349: }
350: return(0);
351: }
353: PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[])
354: {
358: DMGetWorkArray(dm, 6, MPIU_INT, cone);
359: return(0);
360: }
362: PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu)
363: {
364: DM_DA *da = (DM_DA *) dm->data;
365: Vec coordinates;
366: PetscSection section;
367: PetscScalar *coords;
368: PetscReal h[3];
369: PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k;
374: DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);
375: if (dim > 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"The following code only works for dim <= 3");
376: h[0] = (xu - xl)/M;
377: h[1] = (yu - yl)/N;
378: h[2] = (zu - zl)/P;
379: DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);
380: DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);
381: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);
382: PetscSectionSetNumFields(section, 1);
383: PetscSectionSetFieldComponents(section, 0, dim);
384: PetscSectionSetChart(section, vStart, vEnd);
385: for (v = vStart; v < vEnd; ++v) {
386: PetscSectionSetDof(section, v, dim);
387: }
388: PetscSectionSetUp(section);
389: PetscSectionGetStorageSize(section, &size);
390: VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);
391: PetscObjectSetName((PetscObject)coordinates,"coordinates");
392: VecGetArray(coordinates, &coords);
393: for (k = 0; k < nVz; ++k) {
394: PetscInt ind[3], d, off;
396: ind[0] = 0;
397: ind[1] = 0;
398: ind[2] = k + da->zs;
399: for (j = 0; j < nVy; ++j) {
400: ind[1] = j + da->ys;
401: for (i = 0; i < nVx; ++i) {
402: const PetscInt vertex = (k*nVy + j)*nVx + i + vStart;
404: PetscSectionGetOffset(section, vertex, &off);
405: ind[0] = i + da->xs;
406: for (d = 0; d < dim; ++d) {
407: coords[off+d] = h[d]*ind[d];
408: }
409: }
410: }
411: }
412: VecRestoreArray(coordinates, &coords);
413: DMSetCoordinateSection(dm, PETSC_DETERMINE, section);
414: DMSetCoordinatesLocal(dm, coordinates);
415: PetscSectionDestroy(§ion);
416: VecDestroy(&coordinates);
417: return(0);
418: }
420: /* ------------------------------------------------------------------- */
422: /*@C
423: DMDAGetArray - Gets a work array for a DMDA
425: Input Parameter:
426: + da - information about my local patch
427: - ghosted - do you want arrays for the ghosted or nonghosted patch
429: Output Parameters:
430: . vptr - array data structured
432: Note: The vector values are NOT initialized and may have garbage in them, so you may need
433: to zero them.
435: Level: advanced
437: .seealso: DMDARestoreArray()
439: @*/
440: PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr)
441: {
443: PetscInt j,i,xs,ys,xm,ym,zs,zm;
444: char *iarray_start;
445: void **iptr = (void**)vptr;
446: DM_DA *dd = (DM_DA*)da->data;
450: if (ghosted) {
451: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
452: if (dd->arrayghostedin[i]) {
453: *iptr = dd->arrayghostedin[i];
454: iarray_start = (char*)dd->startghostedin[i];
455: dd->arrayghostedin[i] = NULL;
456: dd->startghostedin[i] = NULL;
458: goto done;
459: }
460: }
461: xs = dd->Xs;
462: ys = dd->Ys;
463: zs = dd->Zs;
464: xm = dd->Xe-dd->Xs;
465: ym = dd->Ye-dd->Ys;
466: zm = dd->Ze-dd->Zs;
467: } else {
468: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
469: if (dd->arrayin[i]) {
470: *iptr = dd->arrayin[i];
471: iarray_start = (char*)dd->startin[i];
472: dd->arrayin[i] = NULL;
473: dd->startin[i] = NULL;
475: goto done;
476: }
477: }
478: xs = dd->xs;
479: ys = dd->ys;
480: zs = dd->zs;
481: xm = dd->xe-dd->xs;
482: ym = dd->ye-dd->ys;
483: zm = dd->ze-dd->zs;
484: }
486: switch (da->dim) {
487: case 1: {
488: void *ptr;
490: PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);
492: ptr = (void*)(iarray_start - xs*sizeof(PetscScalar));
493: *iptr = (void*)ptr;
494: break;
495: }
496: case 2: {
497: void **ptr;
499: PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);
501: ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*));
502: for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs);
503: *iptr = (void*)ptr;
504: break;
505: }
506: case 3: {
507: void ***ptr,**bptr;
509: PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);
511: ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*));
512: bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**));
513: for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys);
514: for (i=zs; i<zs+zm; i++) {
515: for (j=ys; j<ys+ym; j++) {
516: ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs);
517: }
518: }
519: *iptr = (void*)ptr;
520: break;
521: }
522: default:
523: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",da->dim);
524: }
526: done:
527: /* add arrays to the checked out list */
528: if (ghosted) {
529: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
530: if (!dd->arrayghostedout[i]) {
531: dd->arrayghostedout[i] = *iptr;
532: dd->startghostedout[i] = iarray_start;
533: break;
534: }
535: }
536: } else {
537: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
538: if (!dd->arrayout[i]) {
539: dd->arrayout[i] = *iptr;
540: dd->startout[i] = iarray_start;
541: break;
542: }
543: }
544: }
545: return(0);
546: }
548: /*@C
549: DMDARestoreArray - Restores an array of derivative types for a DMDA
551: Input Parameter:
552: + da - information about my local patch
553: . ghosted - do you want arrays for the ghosted or nonghosted patch
554: - vptr - array data structured
556: Level: advanced
558: .seealso: DMDAGetArray()
560: @*/
561: PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr)
562: {
563: PetscInt i;
564: void **iptr = (void**)vptr,*iarray_start = 0;
565: DM_DA *dd = (DM_DA*)da->data;
569: if (ghosted) {
570: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
571: if (dd->arrayghostedout[i] == *iptr) {
572: iarray_start = dd->startghostedout[i];
573: dd->arrayghostedout[i] = NULL;
574: dd->startghostedout[i] = NULL;
575: break;
576: }
577: }
578: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
579: if (!dd->arrayghostedin[i]) {
580: dd->arrayghostedin[i] = *iptr;
581: dd->startghostedin[i] = iarray_start;
582: break;
583: }
584: }
585: } else {
586: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
587: if (dd->arrayout[i] == *iptr) {
588: iarray_start = dd->startout[i];
589: dd->arrayout[i] = NULL;
590: dd->startout[i] = NULL;
591: break;
592: }
593: }
594: for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) {
595: if (!dd->arrayin[i]) {
596: dd->arrayin[i] = *iptr;
597: dd->startin[i] = iarray_start;
598: break;
599: }
600: }
601: }
602: return(0);
603: }