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;
17: Xe += Xs;
18: 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: {
39: DM_DA *da = (DM_DA *)dm->data;
40: PetscInt i, xs, xe, Xs, Xe;
41: PetscInt j, ys, ye, Ys, Ye;
42: PetscInt cnt = 0, cell[4], ns = 2;
43: PetscInt c, split[] = {0, 1, 3, 2, 3, 1};
45: if (!da->e) {
46: PetscInt corners[4], nn = 0;
50: switch (da->elementtype) {
51: case DMDA_ELEMENT_Q1:
52: da->nen = 4;
53: break;
54: case DMDA_ELEMENT_P1:
55: da->nen = 3;
56: break;
57: default:
58: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown element type %d", da->elementtype);
59: }
60: nn = da->nen;
62: if (da->elementtype == DMDA_ELEMENT_P1) ns = 2;
63: if (da->elementtype == DMDA_ELEMENT_Q1) ns = 1;
64: DMDAGetCorners(dm, &xs, &ys, NULL, &xe, &ye, NULL);
65: DMDAGetGhostCorners(dm, &Xs, &Ys, NULL, &Xe, &Ye, NULL);
66: xe += xs;
67: Xe += Xs;
68: if (xs != Xs) xs -= 1;
69: ye += ys;
70: Ye += Ys;
71: 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: {
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, 0, 1, 7, 4, 1, 2, 3, 7, 1, 2, 7, 6, 1, 4, 5, 7, 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;
132: Xe += Xs;
133: if (xs != Xs) xs -= 1;
134: ye += ys;
135: Ye += Ys;
136: if (ys != Ys) ys -= 1;
137: ze += zs;
138: Ze += Zs;
139: 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 `DMDA` object
188: Output Parameters:
189: + gx - the x index
190: . gy - the y index
191: - gz - the z index
193: Level: intermediate
195: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`
196: @*/
197: PetscErrorCode DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
198: {
199: PetscInt xs, Xs;
200: PetscInt ys, Ys;
201: PetscInt zs, Zs;
202: PetscBool isda;
208: PetscObjectTypeCompare((PetscObject)da, DMDA, &isda);
210: DMDAGetCorners(da, &xs, &ys, &zs, NULL, NULL, NULL);
211: DMDAGetGhostCorners(da, &Xs, &Ys, &Zs, NULL, NULL, NULL);
212: if (xs != Xs) xs -= 1;
213: if (ys != Ys) ys -= 1;
214: if (zs != Zs) zs -= 1;
215: if (gx) *gx = xs;
216: if (gy) *gy = ys;
217: if (gz) *gz = zs;
218: return 0;
219: }
221: /*@
222: DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by `DMDAGetElements()`
224: Not Collective
226: Input Parameter:
227: . da - the `DMDA` object
229: Output Parameters:
230: + mx - number of local elements in x-direction
231: . my - number of local elements in y-direction
232: - mz - number of local elements in z-direction
234: Level: intermediate
236: Note:
237: It returns the same number of elements, irrespective of the `DMDAElementType`
239: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements`
240: @*/
241: PetscErrorCode DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
242: {
243: PetscInt xs, xe, Xs;
244: PetscInt ys, ye, Ys;
245: PetscInt zs, ze, Zs;
246: PetscInt dim;
247: PetscBool isda;
253: PetscObjectTypeCompare((PetscObject)da, DMDA, &isda);
255: DMDAGetCorners(da, &xs, &ys, &zs, &xe, &ye, &ze);
256: DMDAGetGhostCorners(da, &Xs, &Ys, &Zs, NULL, NULL, NULL);
257: xe += xs;
258: if (xs != Xs) xs -= 1;
259: ye += ys;
260: if (ys != Ys) ys -= 1;
261: ze += zs;
262: if (zs != Zs) zs -= 1;
263: if (mx) *mx = 0;
264: if (my) *my = 0;
265: if (mz) *mz = 0;
266: DMGetDimension(da, &dim);
267: switch (dim) {
268: case 3:
269: if (mz) *mz = ze - zs - 1; /* fall through */
270: case 2:
271: if (my) *my = ye - ys - 1; /* fall through */
272: case 1:
273: if (mx) *mx = xe - xs - 1;
274: break;
275: }
276: return 0;
277: }
279: /*@
280: DMDASetElementType - Sets the element type to be returned by `DMDAGetElements()`
282: Not Collective
284: Input Parameter:
285: . da - the `DMDA` object
287: Output Parameters:
288: . etype - the element type, currently either `DMDA_ELEMENT_P1` or `DMDA_ELEMENT_Q1`
290: Level: intermediate
292: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDAGetElementType()`, `DMDAGetElements()`, `DMDARestoreElements()`
293: @*/
294: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
295: {
296: DM_DA *dd = (DM_DA *)da->data;
297: PetscBool isda;
301: PetscObjectTypeCompare((PetscObject)da, DMDA, &isda);
302: if (!isda) return 0;
303: if (dd->elementtype != etype) {
304: PetscFree(dd->e);
305: ISDestroy(&dd->ecorners);
307: dd->elementtype = etype;
308: dd->ne = 0;
309: dd->nen = 0;
310: dd->e = NULL;
311: }
312: return 0;
313: }
315: /*@
316: DMDAGetElementType - Gets the element type to be returned by `DMDAGetElements()`
318: Not Collective
320: Input Parameter:
321: . da - the `DMDA` object
323: Output Parameters:
324: . etype - the element type, currently either `DMDA_ELEMENT_P1` or `DMDA_ELEMENT_Q1`
326: Level: intermediate
328: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`, `DMDARestoreElements()`
329: @*/
330: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
331: {
332: DM_DA *dd = (DM_DA *)da->data;
333: PetscBool isda;
337: PetscObjectTypeCompare((PetscObject)da, DMDA, &isda);
339: *etype = dd->elementtype;
340: return 0;
341: }
343: /*@C
344: DMDAGetElements - Gets an array containing the indices (in local coordinates)
345: of all the local elements
347: Not Collective
349: Input Parameter:
350: . dm - the `DMDA` object
352: Output Parameters:
353: + nel - number of local elements
354: . nen - number of element nodes
355: - e - the local indices of the elements' vertices
357: Level: intermediate
359: Notes:
360: Call `DMDARestoreElements()` once you have finished accessing the elements.
362: Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.
364: If on each process you integrate over its owned elements and use `ADD_VALUES` in `Vec`/`MatSetValuesLocal()` then you'll obtain the correct result.
366: Fortran Note:
367: Not supported in Fortran
369: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `VecSetValuesLocal()`, `MatSetValuesLocal()`, `DMGlobalToLocalBegin()`, `DMLocalToGlobalBegin()`
370: @*/
371: PetscErrorCode DMDAGetElements(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[])
372: {
373: PetscInt dim;
374: DM_DA *dd = (DM_DA *)dm->data;
375: PetscBool isda;
381: PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda);
384: DMGetDimension(dm, &dim);
385: if (dd->e) {
386: *nel = dd->ne;
387: *nen = dd->nen;
388: *e = dd->e;
389: return 0;
390: }
391: if (dim == -1) {
392: *nel = 0;
393: *nen = 0;
394: *e = NULL;
395: } else if (dim == 1) {
396: DMDAGetElements_1D(dm, nel, nen, e);
397: } else if (dim == 2) {
398: DMDAGetElements_2D(dm, nel, nen, e);
399: } else if (dim == 3) {
400: DMDAGetElements_3D(dm, nel, nen, e);
401: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "DMDA dimension not 1, 2, or 3, it is %" PetscInt_FMT, dim);
402: return 0;
403: }
405: /*@
406: DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates)
407: of the non-overlapping decomposition identified by `DMDAGetElements()`
409: Not Collective
411: Input Parameter:
412: . dm - the `DMDA` object
414: Output Parameters:
415: . is - the index set
417: Level: intermediate
419: Note:
420: Call `DMDARestoreSubdomainCornersIS()` once you have finished accessing the index set.
422: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`, `DMDARestoreElementsCornersIS()`
423: @*/
424: PetscErrorCode DMDAGetSubdomainCornersIS(DM dm, IS *is)
425: {
426: DM_DA *dd = (DM_DA *)dm->data;
427: PetscBool isda;
431: PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda);
434: if (!dd->ecorners) { /* compute elements if not yet done */
435: const PetscInt *e;
436: PetscInt nel, nen;
438: DMDAGetElements(dm, &nel, &nen, &e);
439: DMDARestoreElements(dm, &nel, &nen, &e);
440: }
441: *is = dd->ecorners;
442: return 0;
443: }
445: /*@C
446: DMDARestoreElements - Restores the array obtained with `DMDAGetElements()`
448: Not Collective
450: Input Parameters:
451: + dm - the DM object
452: . nel - number of local elements
453: . nen - number of element nodes
454: - e - the local indices of the elements' vertices
456: Level: intermediate
458: Note:
459: This restore signals the `DMDA` object that you no longer need access to the array information.
461: Fortran Note:
462: Not supported in Fortran
464: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`
465: @*/
466: PetscErrorCode DMDARestoreElements(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[])
467: {
472: *nel = 0;
473: *nen = -1;
474: *e = NULL;
475: return 0;
476: }
478: /*@
479: DMDARestoreSubdomainCornersIS - Restores the `IS` obtained with `DMDAGetSubdomainCornersIS()`
481: Not Collective
483: Input Parameters:
484: + dm - the `DM` object
485: - is - the index set
487: Level: intermediate
489: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetSubdomainCornersIS()`
490: @*/
491: PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm, IS *is)
492: {
495: *is = NULL;
496: return 0;
497: }