Actual source code: dagetelem.c
2: #include <petsc/private/dmdaimpl.h>
4: static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
5: {
6: DM_DA *da = (DM_DA*)dm->data;
7: PetscInt i,xs,xe,Xs,Xe;
8: PetscInt cnt=0;
10: if (!da->e) {
11: PetscInt corners[2];
14: DMDAGetCorners(dm,&xs,NULL,NULL,&xe,NULL,NULL);
15: DMDAGetGhostCorners(dm,&Xs,NULL,NULL,&Xe,NULL,NULL);
16: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
17: da->ne = 1*(xe - xs - 1);
18: PetscMalloc1(1 + 2*da->ne,&da->e);
19: for (i=xs; i<xe-1; i++) {
20: da->e[cnt++] = (i-Xs);
21: da->e[cnt++] = (i-Xs+1);
22: }
23: da->nen = 2;
25: corners[0] = (xs -Xs);
26: corners[1] = (xe-1-Xs);
27: ISCreateGeneral(PETSC_COMM_SELF,2,corners,PETSC_COPY_VALUES,&da->ecorners);
28: }
29: *nel = da->ne;
30: *nen = da->nen;
31: *e = da->e;
32: return 0;
33: }
35: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
36: {
37: DM_DA *da = (DM_DA*)dm->data;
38: PetscInt i,xs,xe,Xs,Xe;
39: PetscInt j,ys,ye,Ys,Ye;
40: PetscInt cnt=0, cell[4], ns=2;
41: PetscInt c, split[] = {0,1,3,
42: 2,3,1};
44: if (!da->e) {
45: PetscInt corners[4],nn = 0;
49: switch (da->elementtype) {
50: case DMDA_ELEMENT_Q1:
51: da->nen = 4;
52: break;
53: case DMDA_ELEMENT_P1:
54: da->nen = 3;
55: break;
56: default:
57: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unknown element type %d",da->elementtype);
58: }
59: nn = da->nen;
61: if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
62: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
63: DMDAGetCorners(dm,&xs,&ys,NULL,&xe,&ye,NULL);
64: DMDAGetGhostCorners(dm,&Xs,&Ys,NULL,&Xe,&Ye,NULL);
65: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
66: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
67: da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
68: PetscMalloc1(1 + nn*da->ne,&da->e);
69: for (j=ys; j<ye-1; j++) {
70: for (i=xs; i<xe-1; i++) {
71: cell[0] = (i-Xs) + (j-Ys)*(Xe-Xs);
72: cell[1] = (i-Xs+1) + (j-Ys)*(Xe-Xs);
73: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
74: cell[3] = (i-Xs) + (j-Ys+1)*(Xe-Xs);
75: if (da->elementtype == DMDA_ELEMENT_P1) {
76: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
77: }
78: if (da->elementtype == DMDA_ELEMENT_Q1) {
79: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
80: }
81: }
82: }
84: corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs);
85: corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs);
86: corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs);
87: corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs);
88: ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);
89: }
90: *nel = da->ne;
91: *nen = da->nen;
92: *e = da->e;
93: return 0;
94: }
96: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
97: {
98: DM_DA *da = (DM_DA*)dm->data;
99: PetscInt i,xs,xe,Xs,Xe;
100: PetscInt j,ys,ye,Ys,Ye;
101: PetscInt k,zs,ze,Zs,Ze;
102: PetscInt cnt=0, cell[8], ns=6;
103: PetscInt c, split[] = {0,1,3,7,
104: 0,1,7,4,
105: 1,2,3,7,
106: 1,2,7,6,
107: 1,4,5,7,
108: 1,5,6,7};
110: if (!da->e) {
111: PetscInt corners[8],nn = 0;
115: switch (da->elementtype) {
116: case DMDA_ELEMENT_Q1:
117: da->nen = 8;
118: break;
119: case DMDA_ELEMENT_P1:
120: da->nen = 4;
121: break;
122: default:
123: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unknown element type %d",da->elementtype);
124: }
125: nn = da->nen;
127: if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;}
128: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
129: DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
130: DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
131: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
132: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
133: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
134: da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
135: PetscMalloc1(1 + nn*da->ne,&da->e);
136: for (k=zs; k<ze-1; k++) {
137: for (j=ys; j<ye-1; j++) {
138: for (i=xs; i<xe-1; i++) {
139: cell[0] = (i -Xs) + (j -Ys)*(Xe-Xs) + (k -Zs)*(Xe-Xs)*(Ye-Ys);
140: cell[1] = (i+1-Xs) + (j -Ys)*(Xe-Xs) + (k -Zs)*(Xe-Xs)*(Ye-Ys);
141: cell[2] = (i+1-Xs) + (j+1-Ys)*(Xe-Xs) + (k -Zs)*(Xe-Xs)*(Ye-Ys);
142: cell[3] = (i -Xs) + (j+1-Ys)*(Xe-Xs) + (k -Zs)*(Xe-Xs)*(Ye-Ys);
143: cell[4] = (i -Xs) + (j -Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
144: cell[5] = (i+1-Xs) + (j -Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
145: cell[6] = (i+1-Xs) + (j+1-Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
146: cell[7] = (i -Xs) + (j+1-Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
147: if (da->elementtype == DMDA_ELEMENT_P1) {
148: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
149: }
150: if (da->elementtype == DMDA_ELEMENT_Q1) {
151: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
152: }
153: }
154: }
155: }
157: corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (zs- Zs)*(Xe-Xs)*(Ye-Ys);
158: corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (zs- Zs)*(Xe-Xs)*(Ye-Ys);
159: corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs- Zs)*(Xe-Xs)*(Ye-Ys);
160: corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs- Zs)*(Xe-Xs)*(Ye-Ys);
161: corners[4] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
162: corners[5] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
163: corners[6] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
164: corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
165: ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);
166: }
167: *nel = da->ne;
168: *nen = da->nen;
169: *e = da->e;
170: return 0;
171: }
173: /*@
174: DMDAGetElementsCorners - Returns the global (x,y,z) indices of the lower left
175: corner of the non-overlapping decomposition identified by DMDAGetElements()
177: Not Collective
179: Input Parameter:
180: . da - the DM object
182: Output Parameters:
183: + gx - the x index
184: . gy - the y index
185: - gz - the z index
187: Level: intermediate
189: Notes:
191: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
192: @*/
193: PetscErrorCode DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
194: {
195: PetscInt xs,Xs;
196: PetscInt ys,Ys;
197: PetscInt zs,Zs;
198: PetscBool isda;
204: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
206: DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);
207: DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
208: if (xs != Xs) xs -= 1;
209: if (ys != Ys) ys -= 1;
210: if (zs != Zs) zs -= 1;
211: if (gx) *gx = xs;
212: if (gy) *gy = ys;
213: if (gz) *gz = zs;
214: return 0;
215: }
217: /*@
218: DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by DMDAGetElements()
220: Not Collective
222: Input Parameter:
223: . da - the DM object
225: Output Parameters:
226: + mx - number of local elements in x-direction
227: . my - number of local elements in y-direction
228: - mz - number of local elements in z-direction
230: Level: intermediate
232: Notes:
233: It returns the same number of elements, irrespective of the DMDAElementType
235: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements
236: @*/
237: PetscErrorCode DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
238: {
239: PetscInt xs,xe,Xs;
240: PetscInt ys,ye,Ys;
241: PetscInt zs,ze,Zs;
242: PetscInt dim;
243: PetscBool isda;
249: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
251: DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);
252: DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
253: xe += xs; if (xs != Xs) xs -= 1;
254: ye += ys; if (ys != Ys) ys -= 1;
255: ze += zs; if (zs != Zs) zs -= 1;
256: if (mx) *mx = 0;
257: if (my) *my = 0;
258: if (mz) *mz = 0;
259: DMGetDimension(da,&dim);
260: switch (dim) {
261: case 3:
262: if (mz) *mz = ze - zs - 1;
263: case 2:
264: if (my) *my = ye - ys - 1;
265: case 1:
266: if (mx) *mx = xe - xs - 1;
267: break;
268: }
269: return 0;
270: }
272: /*@
273: DMDASetElementType - Sets the element type to be returned by DMDAGetElements()
275: Not Collective
277: Input Parameter:
278: . da - the DMDA object
280: Output Parameters:
281: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
283: Level: intermediate
285: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
286: @*/
287: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
288: {
289: DM_DA *dd = (DM_DA*)da->data;
290: PetscBool isda;
294: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
295: if (!isda) return 0;
296: if (dd->elementtype != etype) {
297: PetscFree(dd->e);
298: ISDestroy(&dd->ecorners);
300: dd->elementtype = etype;
301: dd->ne = 0;
302: dd->nen = 0;
303: dd->e = NULL;
304: }
305: return 0;
306: }
308: /*@
309: DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()
311: Not Collective
313: Input Parameter:
314: . da - the DMDA object
316: Output Parameters:
317: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
319: Level: intermediate
321: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
322: @*/
323: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
324: {
325: DM_DA *dd = (DM_DA*)da->data;
326: PetscBool isda;
330: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
332: *etype = dd->elementtype;
333: return 0;
334: }
336: /*@C
337: DMDAGetElements - Gets an array containing the indices (in local coordinates)
338: of all the local elements
340: Not Collective
342: Input Parameter:
343: . dm - the DM object
345: Output Parameters:
346: + nel - number of local elements
347: . nen - number of element nodes
348: - e - the local indices of the elements' vertices
350: Level: intermediate
352: Notes:
353: Call DMDARestoreElements() once you have finished accessing the elements.
355: Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.
357: If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result.
359: Not supported in Fortran
361: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
362: @*/
363: PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
364: {
365: PetscInt dim;
366: DM_DA *dd = (DM_DA*)dm->data;
367: PetscBool isda;
373: PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
376: DMGetDimension(dm, &dim);
377: if (dd->e) {
378: *nel = dd->ne;
379: *nen = dd->nen;
380: *e = dd->e;
381: return 0;
382: }
383: if (dim==-1) {
384: *nel = 0; *nen = 0; *e = NULL;
385: } else if (dim==1) {
386: DMDAGetElements_1D(dm,nel,nen,e);
387: } else if (dim==2) {
388: DMDAGetElements_2D(dm,nel,nen,e);
389: } else if (dim==3) {
390: DMDAGetElements_3D(dm,nel,nen,e);
391: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D",dim);
392: return 0;
393: }
395: /*@
396: DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates)
397: of the non-overlapping decomposition identified by DMDAGetElements
399: Not Collective
401: Input Parameter:
402: . dm - the DM object
404: Output Parameters:
405: . is - the index set
407: Level: intermediate
409: Notes:
410: Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set.
412: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS()
413: @*/
414: PetscErrorCode DMDAGetSubdomainCornersIS(DM dm,IS *is)
415: {
416: DM_DA *dd = (DM_DA*)dm->data;
417: PetscBool isda;
421: PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
424: if (!dd->ecorners) { /* compute elements if not yet done */
425: const PetscInt *e;
426: PetscInt nel,nen;
428: DMDAGetElements(dm,&nel,&nen,&e);
429: DMDARestoreElements(dm,&nel,&nen,&e);
430: }
431: *is = dd->ecorners;
432: return 0;
433: }
435: /*@C
436: DMDARestoreElements - Restores the array obtained with DMDAGetElements()
438: Not Collective
440: Input Parameters:
441: + dm - the DM object
442: . nel - number of local elements
443: . nen - number of element nodes
444: - e - the local indices of the elements' vertices
446: Level: intermediate
448: Note: You should not access these values after you have called this routine.
450: This restore signals the DMDA object that you no longer need access to the array information.
452: Not supported in Fortran
454: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
455: @*/
456: PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
457: {
462: *nel = 0;
463: *nen = -1;
464: *e = NULL;
465: return 0;
466: }
468: /*@
469: DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS()
471: Not Collective
473: Input Parameters:
474: + dm - the DM object
475: - is - the index set
477: Level: intermediate
479: Note:
481: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS()
482: @*/
483: PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm,IS *is)
484: {
487: *is = NULL;
488: return 0;
489: }