Actual source code: dagetelem.c
petsc-3.13.6 2020-09-29
2: #include <petsc/private/dmdaimpl.h>
4: static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
5: {
7: DM_DA *da = (DM_DA*)dm->data;
8: PetscInt i,xs,xe,Xs,Xe;
9: PetscInt cnt=0;
12: if (!da->e) {
13: PetscInt corners[2];
15: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
16: DMDAGetCorners(dm,&xs,0,0,&xe,0,0);
17: DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);
18: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
19: da->ne = 1*(xe - xs - 1);
20: PetscMalloc1(1 + 2*da->ne,&da->e);
21: for (i=xs; i<xe-1; i++) {
22: da->e[cnt++] = (i-Xs);
23: da->e[cnt++] = (i-Xs+1);
24: }
25: da->nen = 2;
27: corners[0] = (xs -Xs);
28: corners[1] = (xe-1-Xs);
29: ISCreateGeneral(PETSC_COMM_SELF,2,corners,PETSC_COPY_VALUES,&da->ecorners);
30: }
31: *nel = da->ne;
32: *nen = da->nen;
33: *e = da->e;
34: return(0);
35: }
37: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
38: {
40: DM_DA *da = (DM_DA*)dm->data;
41: PetscInt i,xs,xe,Xs,Xe;
42: PetscInt j,ys,ye,Ys,Ye;
43: PetscInt cnt=0, cell[4], ns=2;
44: PetscInt c, split[] = {0,1,3,
45: 2,3,1};
48: if (!da->e) {
49: PetscInt corners[4],nn = 0;
51: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
53: switch (da->elementtype) {
54: case DMDA_ELEMENT_Q1:
55: da->nen = 4;
56: break;
57: case DMDA_ELEMENT_P1:
58: da->nen = 3;
59: break;
60: default:
61: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unknown element type %d",da->elementtype);
62: break;
63: }
64: nn = da->nen;
66: if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
67: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
68: DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);
69: DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);
70: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
71: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
72: da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
73: PetscMalloc1(1 + nn*da->ne,&da->e);
74: for (j=ys; j<ye-1; j++) {
75: for (i=xs; i<xe-1; i++) {
76: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs);
77: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs);
78: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
79: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs);
80: if (da->elementtype == DMDA_ELEMENT_P1) {
81: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
82: }
83: if (da->elementtype == DMDA_ELEMENT_Q1) {
84: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
85: }
86: }
87: }
89: corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs);
90: corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs);
91: corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs);
92: corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs);
93: ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);
94: }
95: *nel = da->ne;
96: *nen = da->nen;
97: *e = da->e;
98: return(0);
99: }
101: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
102: {
104: DM_DA *da = (DM_DA*)dm->data;
105: PetscInt i,xs,xe,Xs,Xe;
106: PetscInt j,ys,ye,Ys,Ye;
107: PetscInt k,zs,ze,Zs,Ze;
108: PetscInt cnt=0, cell[8], ns=6;
109: PetscInt c, split[] = {0,1,3,7,
110: 0,1,7,4,
111: 1,2,3,7,
112: 1,2,7,6,
113: 1,4,5,7,
114: 1,5,6,7};
117: if (!da->e) {
118: PetscInt corners[8],nn = 0;
120: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
122: switch (da->elementtype) {
123: case DMDA_ELEMENT_Q1:
124: da->nen = 8;
125: break;
126: case DMDA_ELEMENT_P1:
127: da->nen = 4;
128: break;
129: default:
130: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unknown element type %d",da->elementtype);
131: break;
132: }
133: nn = da->nen;
135: if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;}
136: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
137: DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
138: DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
139: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
140: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
141: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
142: da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
143: PetscMalloc1(1 + nn*da->ne,&da->e);
144: for (k=zs; k<ze-1; k++) {
145: for (j=ys; j<ye-1; j++) {
146: for (i=xs; i<xe-1; i++) {
147: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
148: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
149: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
150: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
151: cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
152: cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
153: cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
154: cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
155: if (da->elementtype == DMDA_ELEMENT_P1) {
156: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
157: }
158: if (da->elementtype == DMDA_ELEMENT_Q1) {
159: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
160: }
161: }
162: }
163: }
165: corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys);
166: corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys);
167: corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys);
168: corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys);
169: corners[4] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
170: corners[5] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
171: corners[6] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
172: corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
173: ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);
174: }
175: *nel = da->ne;
176: *nen = da->nen;
177: *e = da->e;
178: return(0);
179: }
181: /*@
182: DMDAGetElementsCorners - Returns the global (x,y,z) indices of the lower left
183: corner of the non-overlapping decomposition identified by DMDAGetElements()
185: Not Collective
187: Input Parameter:
188: . da - the DM object
190: Output Parameters:
191: + gx - the x index
192: . gy - the y index
193: - gz - the z index
195: Level: intermediate
197: Notes:
199: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
200: @*/
201: PetscErrorCode DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
202: {
203: PetscInt xs,Xs;
204: PetscInt ys,Ys;
205: PetscInt zs,Zs;
206: PetscBool isda;
214: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
215: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
216: DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);
217: DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
218: if (xs != Xs) xs -= 1;
219: if (ys != Ys) ys -= 1;
220: if (zs != Zs) zs -= 1;
221: if (gx) *gx = xs;
222: if (gy) *gy = ys;
223: if (gz) *gz = zs;
224: return(0);
225: }
227: /*@
228: DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by DMDAGetElements()
230: Not Collective
232: Input Parameter:
233: . da - the DM object
235: Output Parameters:
236: + mx - number of local elements in x-direction
237: . my - number of local elements in y-direction
238: - mz - number of local elements in z-direction
240: Level: intermediate
242: Notes:
243: It returns the same number of elements, irrespective of the DMDAElementType
245: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements
246: @*/
247: PetscErrorCode DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
248: {
249: PetscInt xs,xe,Xs;
250: PetscInt ys,ye,Ys;
251: PetscInt zs,ze,Zs;
252: PetscInt dim;
253: PetscBool isda;
261: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
262: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
263: DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);
264: DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
265: xe += xs; if (xs != Xs) xs -= 1;
266: ye += ys; if (ys != Ys) ys -= 1;
267: ze += zs; if (zs != Zs) zs -= 1;
268: if (mx) *mx = 0;
269: if (my) *my = 0;
270: if (mz) *mz = 0;
271: DMGetDimension(da,&dim);
272: switch (dim) {
273: case 3:
274: if (mz) *mz = ze - zs - 1;
275: case 2:
276: if (my) *my = ye - ys - 1;
277: case 1:
278: if (mx) *mx = xe - xs - 1;
279: break;
280: }
281: return(0);
282: }
284: /*@
285: DMDASetElementType - Sets the element type to be returned by DMDAGetElements()
287: Not Collective
289: Input Parameter:
290: . da - the DMDA object
292: Output Parameters:
293: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
295: Level: intermediate
297: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
298: @*/
299: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
300: {
301: DM_DA *dd = (DM_DA*)da->data;
303: PetscBool isda;
308: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
309: if (!isda) return(0);
310: if (dd->elementtype != etype) {
311: PetscFree(dd->e);
312: ISDestroy(&dd->ecorners);
314: dd->elementtype = etype;
315: dd->ne = 0;
316: dd->nen = 0;
317: dd->e = NULL;
318: }
319: return(0);
320: }
322: /*@
323: DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()
325: Not Collective
327: Input Parameter:
328: . da - the DMDA object
330: Output Parameters:
331: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
333: Level: intermediate
335: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
336: @*/
337: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
338: {
339: DM_DA *dd = (DM_DA*)da->data;
341: PetscBool isda;
346: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
347: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
348: *etype = dd->elementtype;
349: return(0);
350: }
352: /*@C
353: DMDAGetElements - Gets an array containing the indices (in local coordinates)
354: of all the local elements
356: Not Collective
358: Input Parameter:
359: . dm - the DM object
361: Output Parameters:
362: + nel - number of local elements
363: . nen - number of element nodes
364: - e - the local indices of the elements' vertices
366: Level: intermediate
368: Notes:
369: Call DMDARestoreElements() once you have finished accessing the elements.
371: Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.
373: If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result.
375: Not supported in Fortran
377: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
378: @*/
379: PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
380: {
381: PetscInt dim;
383: DM_DA *dd = (DM_DA*)dm->data;
384: PetscBool isda;
391: PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
392: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
393: if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElements() requires you use a stencil type of DMDA_STENCIL_BOX");
394: DMGetDimension(dm, &dim);
395: if (dd->e) {
396: *nel = dd->ne;
397: *nen = dd->nen;
398: *e = dd->e;
399: return(0);
400: }
401: if (dim==-1) {
402: *nel = 0; *nen = 0; *e = NULL;
403: } else if (dim==1) {
404: DMDAGetElements_1D(dm,nel,nen,e);
405: } else if (dim==2) {
406: DMDAGetElements_2D(dm,nel,nen,e);
407: } else if (dim==3) {
408: DMDAGetElements_3D(dm,nel,nen,e);
409: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
410: return(0);
411: }
413: /*@
414: DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates)
415: of the non-overlapping decomposition identified by DMDAGetElements
417: Not Collective
419: Input Parameter:
420: . dm - the DM object
422: Output Parameters:
423: . is - the index set
425: Level: intermediate
427: Notes:
428: Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set.
430: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS()
431: @*/
432: PetscErrorCode DMDAGetSubdomainCornersIS(DM dm,IS *is)
433: {
435: DM_DA *dd = (DM_DA*)dm->data;
436: PetscBool isda;
441: PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
442: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
443: if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElement() requires you use a stencil type of DMDA_STENCIL_BOX");
444: if (!dd->ecorners) { /* compute elements if not yet done */
445: const PetscInt *e;
446: PetscInt nel,nen;
448: DMDAGetElements(dm,&nel,&nen,&e);
449: DMDARestoreElements(dm,&nel,&nen,&e);
450: }
451: *is = dd->ecorners;
452: return(0);
453: }
455: /*@C
456: DMDARestoreElements - Restores the array obtained with DMDAGetElements()
458: Not Collective
460: Input Parameter:
461: + dm - the DM object
462: . nel - number of local elements
463: . nen - number of element nodes
464: - e - the local indices of the elements' vertices
466: Level: intermediate
468: Note: You should not access these values after you have called this routine.
470: This restore signals the DMDA object that you no longer need access to the array information.
472: Not supported in Fortran
474: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
475: @*/
476: PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
477: {
483: *nel = 0;
484: *nen = -1;
485: *e = NULL;
486: return(0);
487: }
489: /*@
490: DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS()
492: Not Collective
494: Input Parameter:
495: + dm - the DM object
496: - is - the index set
498: Level: intermediate
500: Note:
502: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS()
503: @*/
504: PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm,IS *is)
505: {
509: *is = NULL;
510: return(0);
511: }