Actual source code: dagetelem.c
petsc-3.14.6 2021-03-30
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,NULL,NULL,&xe,NULL,NULL);
17: DMDAGetGhostCorners(dm,&Xs,NULL,NULL,&Xe,NULL,NULL);
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: }
63: nn = da->nen;
65: if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
66: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
67: DMDAGetCorners(dm,&xs,&ys,NULL,&xe,&ye,NULL);
68: DMDAGetGhostCorners(dm,&Xs,&Ys,NULL,&Xe,&Ye,NULL);
69: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
70: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
71: da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
72: PetscMalloc1(1 + nn*da->ne,&da->e);
73: for (j=ys; j<ye-1; j++) {
74: for (i=xs; i<xe-1; i++) {
75: cell[0] = (i-Xs) + (j-Ys)*(Xe-Xs);
76: cell[1] = (i-Xs+1) + (j-Ys)*(Xe-Xs);
77: cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
78: cell[3] = (i-Xs) + (j-Ys+1)*(Xe-Xs);
79: if (da->elementtype == DMDA_ELEMENT_P1) {
80: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
81: }
82: if (da->elementtype == DMDA_ELEMENT_Q1) {
83: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
84: }
85: }
86: }
88: corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs);
89: corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs);
90: corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs);
91: corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs);
92: ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);
93: }
94: *nel = da->ne;
95: *nen = da->nen;
96: *e = da->e;
97: return(0);
98: }
100: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
101: {
103: DM_DA *da = (DM_DA*)dm->data;
104: PetscInt i,xs,xe,Xs,Xe;
105: PetscInt j,ys,ye,Ys,Ye;
106: PetscInt k,zs,ze,Zs,Ze;
107: PetscInt cnt=0, cell[8], ns=6;
108: PetscInt c, split[] = {0,1,3,7,
109: 0,1,7,4,
110: 1,2,3,7,
111: 1,2,7,6,
112: 1,4,5,7,
113: 1,5,6,7};
116: if (!da->e) {
117: PetscInt corners[8],nn = 0;
119: if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
121: switch (da->elementtype) {
122: case DMDA_ELEMENT_Q1:
123: da->nen = 8;
124: break;
125: case DMDA_ELEMENT_P1:
126: da->nen = 4;
127: break;
128: default:
129: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unknown element type %d",da->elementtype);
130: }
131: nn = da->nen;
133: if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;}
134: if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
135: DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
136: DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
137: xe += xs; Xe += Xs; if (xs != Xs) xs -= 1;
138: ye += ys; Ye += Ys; if (ys != Ys) ys -= 1;
139: ze += zs; Ze += Zs; if (zs != Zs) zs -= 1;
140: da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
141: PetscMalloc1(1 + nn*da->ne,&da->e);
142: for (k=zs; k<ze-1; k++) {
143: for (j=ys; j<ye-1; j++) {
144: for (i=xs; i<xe-1; i++) {
145: cell[0] = (i -Xs) + (j -Ys)*(Xe-Xs) + (k -Zs)*(Xe-Xs)*(Ye-Ys);
146: cell[1] = (i+1-Xs) + (j -Ys)*(Xe-Xs) + (k -Zs)*(Xe-Xs)*(Ye-Ys);
147: cell[2] = (i+1-Xs) + (j+1-Ys)*(Xe-Xs) + (k -Zs)*(Xe-Xs)*(Ye-Ys);
148: cell[3] = (i -Xs) + (j+1-Ys)*(Xe-Xs) + (k -Zs)*(Xe-Xs)*(Ye-Ys);
149: cell[4] = (i -Xs) + (j -Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
150: cell[5] = (i+1-Xs) + (j -Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
151: cell[6] = (i+1-Xs) + (j+1-Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
152: cell[7] = (i -Xs) + (j+1-Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
153: if (da->elementtype == DMDA_ELEMENT_P1) {
154: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
155: }
156: if (da->elementtype == DMDA_ELEMENT_Q1) {
157: for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
158: }
159: }
160: }
161: }
163: corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (zs- Zs)*(Xe-Xs)*(Ye-Ys);
164: corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (zs- Zs)*(Xe-Xs)*(Ye-Ys);
165: corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs- Zs)*(Xe-Xs)*(Ye-Ys);
166: corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs- Zs)*(Xe-Xs)*(Ye-Ys);
167: corners[4] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
168: corners[5] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
169: corners[6] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
170: corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
171: ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);
172: }
173: *nel = da->ne;
174: *nen = da->nen;
175: *e = da->e;
176: return(0);
177: }
179: /*@
180: DMDAGetElementsCorners - Returns the global (x,y,z) indices of the lower left
181: corner of the non-overlapping decomposition identified by DMDAGetElements()
183: Not Collective
185: Input Parameter:
186: . da - the DM object
188: Output Parameters:
189: + gx - the x index
190: . gy - the y index
191: - gz - the z index
193: Level: intermediate
195: Notes:
197: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
198: @*/
199: PetscErrorCode DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
200: {
201: PetscInt xs,Xs;
202: PetscInt ys,Ys;
203: PetscInt zs,Zs;
204: PetscBool isda;
212: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
213: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
214: DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);
215: DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
216: if (xs != Xs) xs -= 1;
217: if (ys != Ys) ys -= 1;
218: if (zs != Zs) zs -= 1;
219: if (gx) *gx = xs;
220: if (gy) *gy = ys;
221: if (gz) *gz = zs;
222: return(0);
223: }
225: /*@
226: DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by DMDAGetElements()
228: Not Collective
230: Input Parameter:
231: . da - the DM object
233: Output Parameters:
234: + mx - number of local elements in x-direction
235: . my - number of local elements in y-direction
236: - mz - number of local elements in z-direction
238: Level: intermediate
240: Notes:
241: It returns the same number of elements, irrespective of the DMDAElementType
243: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements
244: @*/
245: PetscErrorCode DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
246: {
247: PetscInt xs,xe,Xs;
248: PetscInt ys,ye,Ys;
249: PetscInt zs,ze,Zs;
250: PetscInt dim;
251: PetscBool isda;
259: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
260: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
261: DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);
262: DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
263: xe += xs; if (xs != Xs) xs -= 1;
264: ye += ys; if (ys != Ys) ys -= 1;
265: ze += zs; if (zs != Zs) zs -= 1;
266: if (mx) *mx = 0;
267: if (my) *my = 0;
268: if (mz) *mz = 0;
269: DMGetDimension(da,&dim);
270: switch (dim) {
271: case 3:
272: if (mz) *mz = ze - zs - 1;
273: case 2:
274: if (my) *my = ye - ys - 1;
275: case 1:
276: if (mx) *mx = xe - xs - 1;
277: break;
278: }
279: return(0);
280: }
282: /*@
283: DMDASetElementType - Sets the element type to be returned by DMDAGetElements()
285: Not Collective
287: Input Parameter:
288: . da - the DMDA object
290: Output Parameters:
291: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
293: Level: intermediate
295: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
296: @*/
297: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
298: {
299: DM_DA *dd = (DM_DA*)da->data;
301: PetscBool isda;
306: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
307: if (!isda) return(0);
308: if (dd->elementtype != etype) {
309: PetscFree(dd->e);
310: ISDestroy(&dd->ecorners);
312: dd->elementtype = etype;
313: dd->ne = 0;
314: dd->nen = 0;
315: dd->e = NULL;
316: }
317: return(0);
318: }
320: /*@
321: DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()
323: Not Collective
325: Input Parameter:
326: . da - the DMDA object
328: Output Parameters:
329: . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1
331: Level: intermediate
333: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
334: @*/
335: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
336: {
337: DM_DA *dd = (DM_DA*)da->data;
339: PetscBool isda;
344: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
345: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
346: *etype = dd->elementtype;
347: return(0);
348: }
350: /*@C
351: DMDAGetElements - Gets an array containing the indices (in local coordinates)
352: of all the local elements
354: Not Collective
356: Input Parameter:
357: . dm - the DM object
359: Output Parameters:
360: + nel - number of local elements
361: . nen - number of element nodes
362: - e - the local indices of the elements' vertices
364: Level: intermediate
366: Notes:
367: Call DMDARestoreElements() once you have finished accessing the elements.
369: Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.
371: If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result.
373: Not supported in Fortran
375: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
376: @*/
377: PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
378: {
379: PetscInt dim;
381: DM_DA *dd = (DM_DA*)dm->data;
382: PetscBool isda;
389: PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
390: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
391: 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");
392: DMGetDimension(dm, &dim);
393: if (dd->e) {
394: *nel = dd->ne;
395: *nen = dd->nen;
396: *e = dd->e;
397: return(0);
398: }
399: if (dim==-1) {
400: *nel = 0; *nen = 0; *e = NULL;
401: } else if (dim==1) {
402: DMDAGetElements_1D(dm,nel,nen,e);
403: } else if (dim==2) {
404: DMDAGetElements_2D(dm,nel,nen,e);
405: } else if (dim==3) {
406: DMDAGetElements_3D(dm,nel,nen,e);
407: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
408: return(0);
409: }
411: /*@
412: DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates)
413: of the non-overlapping decomposition identified by DMDAGetElements
415: Not Collective
417: Input Parameter:
418: . dm - the DM object
420: Output Parameters:
421: . is - the index set
423: Level: intermediate
425: Notes:
426: Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set.
428: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS()
429: @*/
430: PetscErrorCode DMDAGetSubdomainCornersIS(DM dm,IS *is)
431: {
433: DM_DA *dd = (DM_DA*)dm->data;
434: PetscBool isda;
439: PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
440: if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
441: 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");
442: if (!dd->ecorners) { /* compute elements if not yet done */
443: const PetscInt *e;
444: PetscInt nel,nen;
446: DMDAGetElements(dm,&nel,&nen,&e);
447: DMDARestoreElements(dm,&nel,&nen,&e);
448: }
449: *is = dd->ecorners;
450: return(0);
451: }
453: /*@C
454: DMDARestoreElements - Restores the array obtained with DMDAGetElements()
456: Not Collective
458: Input Parameter:
459: + dm - the DM object
460: . nel - number of local elements
461: . nen - number of element nodes
462: - e - the local indices of the elements' vertices
464: Level: intermediate
466: Note: You should not access these values after you have called this routine.
468: This restore signals the DMDA object that you no longer need access to the array information.
470: Not supported in Fortran
472: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
473: @*/
474: PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
475: {
481: *nel = 0;
482: *nen = -1;
483: *e = NULL;
484: return(0);
485: }
487: /*@
488: DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS()
490: Not Collective
492: Input Parameter:
493: + dm - the DM object
494: - is - the index set
496: Level: intermediate
498: Note:
500: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS()
501: @*/
502: PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm,IS *is)
503: {
507: *is = NULL;
508: return(0);
509: }