Actual source code: dagetelem.c
petsc-3.9.4 2018-09-11
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: 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 = 2;
31: *e = da->e;
32: return(0);
33: }
35: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
36: {
38: DM_DA *da = (DM_DA*)dm->data;
39: PetscInt i,xs,xe,Xs,Xe;
40: PetscInt j,ys,ye,Ys,Ye;
41: PetscInt cnt=0, cell[4], ns=2, nn=3;
42: PetscInt c, split[] = {0,1,3,
43: 2,3,1};
46: if (da->elementtype == DMDA_ELEMENT_P1) {nn=3;}
47: if (da->elementtype == DMDA_ELEMENT_Q1) {nn=4;}
48: if (!da->e) {
49: PetscInt corners[4];
51: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
52: if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
53: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
54: DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);
55: DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);
56: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
57: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
58: da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
59: PetscMalloc1(1 + nn*da->ne,&da->e);
60: for (j=ys; j<ye-1; j++) {
61: for (i=xs; i<xe-1; i++) {
62: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs);
63: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs);
64: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
65: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs);
66: if (da->elementtype == DMDA_ELEMENT_P1) {
67: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
68: }
69: if (da->elementtype == DMDA_ELEMENT_Q1) {
70: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
71: }
72: }
73: }
74: corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs);
75: corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs);
76: corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs);
77: corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs);
78: ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);
79: }
80: *nel = da->ne;
81: *nen = nn;
82: *e = da->e;
83: return(0);
84: }
86: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
87: {
89: DM_DA *da = (DM_DA*)dm->data;
90: PetscInt i,xs,xe,Xs,Xe;
91: PetscInt j,ys,ye,Ys,Ye;
92: PetscInt k,zs,ze,Zs,Ze;
93: PetscInt cnt=0, cell[8], ns=6, nn=4;
94: PetscInt c, split[] = {0,1,3,7,
95: 0,1,7,4,
96: 1,2,3,7,
97: 1,2,7,6,
98: 1,4,5,7,
99: 1,5,6,7};
102: if (da->elementtype == DMDA_ELEMENT_P1) {nn=4;}
103: if (da->elementtype == DMDA_ELEMENT_Q1) {nn=8;}
104: if (!da->e) {
105: PetscInt corners[8];
107: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
108: if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;}
109: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
110: DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
111: DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
112: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
113: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
114: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
115: da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
116: PetscMalloc1(1 + nn*da->ne,&da->e);
117: for (k=zs; k<ze-1; k++) {
118: for (j=ys; j<ye-1; j++) {
119: for (i=xs; i<xe-1; i++) {
120: cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
121: cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
122: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
123: cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys);
124: cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
125: cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
126: cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
127: cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
128: if (da->elementtype == DMDA_ELEMENT_P1) {
129: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
130: }
131: if (da->elementtype == DMDA_ELEMENT_Q1) {
132: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
133: }
134: }
135: }
136: }
137: corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys);
138: corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys);
139: corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys);
140: corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys);
141: corners[4] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
142: corners[5] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
143: corners[6] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
144: corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
145: ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);
146: }
147: *nel = da->ne;
148: *nen = nn;
149: *e = da->e;
150: return(0);
151: }
153: /*@
154: DMDAGetElementsCorners - Returns the global (x,y,z) indices of the lower left
155: corner of the non-overlapping decomposition identified by DMDAGetElements()
157: Not Collective
159: Input Parameter:
160: . da - the DM object
162: Output Parameters:
163: + gx - the x index
164: . gy - the y index
165: - gz - the z index
167: Level: intermediate
169: Notes:
171: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
172: @*/
173: PetscErrorCode DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
174: {
175: PetscInt xs,Xs;
176: PetscInt ys,Ys;
177: PetscInt zs,Zs;
178: PetscBool isda;
186: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
187: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
188: DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);
189: DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
190: if (xs != Xs) xs -= 1;
191: if (ys != Ys) ys -= 1;
192: if (zs != Zs) zs -= 1;
193: if (gx) *gx = xs;
194: if (gy) *gy = ys;
195: if (gz) *gz = zs;
196: return(0);
197: }
199: /*@
200: DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by DMDAGetElements()
202: Not Collective
204: Input Parameter:
205: . da - the DM object
207: Output Parameters:
208: + mx - number of local elements in x-direction
209: . my - number of local elements in y-direction
210: - mz - number of local elements in z-direction
212: Level: intermediate
214: Notes: It returns the same number of elements, irrespective of the DMDAElementType
216: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements
217: @*/
218: PetscErrorCode DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
219: {
220: PetscInt xs,xe,Xs;
221: PetscInt ys,ye,Ys;
222: PetscInt zs,ze,Zs;
223: PetscInt dim;
224: PetscBool isda;
232: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
233: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
234: DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);
235: DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
236: xe += xs; if (xs != Xs) xs -= 1;
237: ye += ys; if (ys != Ys) ys -= 1;
238: ze += zs; if (zs != Zs) zs -= 1;
239: if (mx) *mx = 0;
240: if (my) *my = 0;
241: if (mz) *mz = 0;
242: DMGetDimension(da,&dim);
243: switch (dim) {
244: case 3:
245: if (mz) *mz = ze - zs - 1;
246: case 2:
247: if (my) *my = ye - ys - 1;
248: case 1:
249: if (mx) *mx = xe - xs - 1;
250: break;
251: }
252: return(0);
253: }
255: /*@
256: DMDASetElementType - Sets the element type to be returned by DMDAGetElements()
258: Not Collective
260: Input Parameter:
261: . da - the DMDA object
263: Output Parameters:
264: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
266: Level: intermediate
268: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
269: @*/
270: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
271: {
272: DM_DA *dd = (DM_DA*)da->data;
274: PetscBool isda;
279: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
280: if (!isda) return(0);
281: if (dd->elementtype != etype) {
282: PetscFree(dd->e);
283: ISDestroy(&dd->ecorners);
285: dd->elementtype = etype;
286: dd->ne = 0;
287: dd->e = NULL;
288: }
289: return(0);
290: }
292: /*@
293: DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()
295: Not Collective
297: Input Parameter:
298: . da - the DMDA object
300: Output Parameters:
301: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
303: Level: intermediate
305: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
306: @*/
307: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
308: {
309: DM_DA *dd = (DM_DA*)da->data;
311: PetscBool isda;
316: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
317: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
318: *etype = dd->elementtype;
319: return(0);
320: }
322: /*@C
323: DMDAGetElements - Gets an array containing the indices (in local coordinates)
324: of all the local elements
326: Not Collective
328: Input Parameter:
329: . dm - the DM object
331: Output Parameters:
332: + nel - number of local elements
333: . nen - number of element nodes
334: - e - the local indices of the elements' vertices
336: Level: intermediate
338: Notes:
339: Call DMDARestoreElements() once you have finished accessing the elements.
341: Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.
343: If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result.
345: Not supported in Fortran
347: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
348: @*/
349: PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
350: {
351: PetscInt dim;
353: DM_DA *dd = (DM_DA*)dm->data;
354: PetscBool isda;
361: PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
362: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
363: 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");
364: DMGetDimension(dm, &dim);
365: if (dim==-1) {
366: *nel = 0; *nen = 0; *e = NULL;
367: } else if (dim==1) {
368: DMDAGetElements_1D(dm,nel,nen,e);
369: } else if (dim==2) {
370: DMDAGetElements_2D(dm,nel,nen,e);
371: } else if (dim==3) {
372: DMDAGetElements_3D(dm,nel,nen,e);
373: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
374: return(0);
375: }
377: /*@
378: DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates)
379: of the non-overlapping decomposition identified by DMDAGetElements
381: Not Collective
383: Input Parameter:
384: . dm - the DM object
386: Output Parameters:
387: . is - the index set
389: Level: intermediate
391: Notes: Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set.
393: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS()
394: @*/
395: PetscErrorCode DMDAGetSubdomainCornersIS(DM dm,IS *is)
396: {
398: DM_DA *dd = (DM_DA*)dm->data;
399: PetscBool isda;
404: PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
405: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
406: 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");
407: if (!dd->ecorners) { /* compute elements if not yet done */
408: const PetscInt *e;
409: PetscInt nel,nen;
411: DMDAGetElements(dm,&nel,&nen,&e);
412: DMDARestoreElements(dm,&nel,&nen,&e);
413: }
414: *is = dd->ecorners;
415: return(0);
416: }
418: /*@C
419: DMDARestoreElements - Restores the array obtained with DMDAGetElements()
421: Not Collective
423: Input Parameter:
424: + dm - the DM object
425: . nel - number of local elements
426: . nen - number of element nodes
427: - e - the local indices of the elements' vertices
429: Level: intermediate
431: Note: You should not access these values after you have called this routine.
433: This restore signals the DMDA object that you no longer need access to the array information.
435: Not supported in Fortran
437: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
438: @*/
439: PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
440: {
446: *nel = 0;
447: *nen = -1;
448: *e = NULL;
449: return(0);
450: }
452: /*@
453: DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS()
455: Not Collective
457: Input Parameter:
458: + dm - the DM object
459: - is - the index set
461: Level: intermediate
463: Note:
465: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS()
466: @*/
467: PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm,IS *is)
468: {
472: *is = NULL;
473: return(0);
474: }