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:   PetscFunctionBegin;
 11:   if (!da->e) {
 12:     PetscInt corners[2];

 14:     PetscCheck(da->s, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot get elements for DMDA with zero stencil width");
 15:     PetscCall(DMDAGetCorners(dm, &xs, NULL, NULL, &xe, NULL, NULL));
 16:     PetscCall(DMDAGetGhostCorners(dm, &Xs, NULL, NULL, &Xe, NULL, NULL));
 17:     xe += xs;
 18:     Xe += Xs;
 19:     if (xs != Xs) xs -= 1;
 20:     da->ne = 1 * (xe - xs - 1);
 21:     PetscCall(PetscMalloc1(1 + 2 * da->ne, &da->e));
 22:     for (i = xs; i < xe - 1; i++) {
 23:       da->e[cnt++] = (i - Xs);
 24:       da->e[cnt++] = (i - Xs + 1);
 25:     }
 26:     da->nen = 2;

 28:     corners[0] = (xs - Xs);
 29:     corners[1] = (xe - 1 - Xs);
 30:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 2, corners, PETSC_COPY_VALUES, &da->ecorners));
 31:   }
 32:   *nel = da->ne;
 33:   *nen = da->nen;
 34:   *e   = da->e;
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }

 38: static PetscErrorCode DMDAGetElements_2D(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[])
 39: {
 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, 2, 3, 1};

 46:   PetscFunctionBegin;
 47:   if (!da->e) {
 48:     PetscInt corners[4], nn = 0;

 50:     PetscCheck(da->s, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot get elements for DMDA with zero stencil width");

 52:     switch (da->elementtype) {
 53:     case DMDA_ELEMENT_Q1:
 54:       da->nen = 4;
 55:       break;
 56:     case DMDA_ELEMENT_P1:
 57:       da->nen = 3;
 58:       break;
 59:     default:
 60:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown element type %d", da->elementtype);
 61:     }
 62:     nn = da->nen;

 64:     if (da->elementtype == DMDA_ELEMENT_P1) ns = 2;
 65:     if (da->elementtype == DMDA_ELEMENT_Q1) ns = 1;
 66:     PetscCall(DMDAGetCorners(dm, &xs, &ys, NULL, &xe, &ye, NULL));
 67:     PetscCall(DMDAGetGhostCorners(dm, &Xs, &Ys, NULL, &Xe, &Ye, NULL));
 68:     xe += xs;
 69:     Xe += Xs;
 70:     if (xs != Xs) xs -= 1;
 71:     ye += ys;
 72:     Ye += Ys;
 73:     if (ys != Ys) ys -= 1;
 74:     da->ne = ns * (xe - xs - 1) * (ye - ys - 1);
 75:     PetscCall(PetscMalloc1(1 + nn * da->ne, &da->e));
 76:     for (j = ys; j < ye - 1; j++) {
 77:       for (i = xs; i < xe - 1; i++) {
 78:         cell[0] = (i - Xs) + (j - Ys) * (Xe - Xs);
 79:         cell[1] = (i - Xs + 1) + (j - Ys) * (Xe - Xs);
 80:         cell[2] = (i - Xs + 1) + (j - Ys + 1) * (Xe - Xs);
 81:         cell[3] = (i - Xs) + (j - Ys + 1) * (Xe - Xs);
 82:         if (da->elementtype == DMDA_ELEMENT_P1) {
 83:           for (c = 0; c < ns * nn; c++) da->e[cnt++] = cell[split[c]];
 84:         }
 85:         if (da->elementtype == DMDA_ELEMENT_Q1) {
 86:           for (c = 0; c < ns * nn; c++) da->e[cnt++] = cell[c];
 87:         }
 88:       }
 89:     }

 91:     corners[0] = (xs - Xs) + (ys - Ys) * (Xe - Xs);
 92:     corners[1] = (xe - 1 - Xs) + (ys - Ys) * (Xe - Xs);
 93:     corners[2] = (xs - Xs) + (ye - 1 - Ys) * (Xe - Xs);
 94:     corners[3] = (xe - 1 - Xs) + (ye - 1 - Ys) * (Xe - Xs);
 95:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 4, corners, PETSC_COPY_VALUES, &da->ecorners));
 96:   }
 97:   *nel = da->ne;
 98:   *nen = da->nen;
 99:   *e   = da->e;
100:   PetscFunctionReturn(PETSC_SUCCESS);
101: }

103: static PetscErrorCode DMDAGetElements_3D(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[])
104: {
105:   DM_DA   *da = (DM_DA *)dm->data;
106:   PetscInt i, xs, xe, Xs, Xe;
107:   PetscInt j, ys, ye, Ys, Ye;
108:   PetscInt k, zs, ze, Zs, Ze;
109:   PetscInt cnt = 0, cell[8], ns = 6;
110:   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};

112:   PetscFunctionBegin;
113:   if (!da->e) {
114:     PetscInt corners[8], nn = 0;

116:     PetscCheck(da->s, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot get elements for DMDA with zero stencil width");

118:     switch (da->elementtype) {
119:     case DMDA_ELEMENT_Q1:
120:       da->nen = 8;
121:       break;
122:     case DMDA_ELEMENT_P1:
123:       da->nen = 4;
124:       break;
125:     default:
126:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unknown element type %d", da->elementtype);
127:     }
128:     nn = da->nen;

130:     if (da->elementtype == DMDA_ELEMENT_P1) ns = 6;
131:     if (da->elementtype == DMDA_ELEMENT_Q1) ns = 1;
132:     PetscCall(DMDAGetCorners(dm, &xs, &ys, &zs, &xe, &ye, &ze));
133:     PetscCall(DMDAGetGhostCorners(dm, &Xs, &Ys, &Zs, &Xe, &Ye, &Ze));
134:     xe += xs;
135:     Xe += Xs;
136:     if (xs != Xs) xs -= 1;
137:     ye += ys;
138:     Ye += Ys;
139:     if (ys != Ys) ys -= 1;
140:     ze += zs;
141:     Ze += Zs;
142:     if (zs != Zs) zs -= 1;
143:     da->ne = ns * (xe - xs - 1) * (ye - ys - 1) * (ze - zs - 1);
144:     PetscCall(PetscMalloc1(1 + nn * da->ne, &da->e));
145:     for (k = zs; k < ze - 1; k++) {
146:       for (j = ys; j < ye - 1; j++) {
147:         for (i = xs; i < xe - 1; i++) {
148:           cell[0] = (i - Xs) + (j - Ys) * (Xe - Xs) + (k - Zs) * (Xe - Xs) * (Ye - Ys);
149:           cell[1] = (i + 1 - Xs) + (j - Ys) * (Xe - Xs) + (k - Zs) * (Xe - Xs) * (Ye - Ys);
150:           cell[2] = (i + 1 - Xs) + (j + 1 - Ys) * (Xe - Xs) + (k - Zs) * (Xe - Xs) * (Ye - Ys);
151:           cell[3] = (i - Xs) + (j + 1 - Ys) * (Xe - Xs) + (k - Zs) * (Xe - Xs) * (Ye - Ys);
152:           cell[4] = (i - Xs) + (j - Ys) * (Xe - Xs) + (k + 1 - Zs) * (Xe - Xs) * (Ye - Ys);
153:           cell[5] = (i + 1 - Xs) + (j - Ys) * (Xe - Xs) + (k + 1 - Zs) * (Xe - Xs) * (Ye - Ys);
154:           cell[6] = (i + 1 - Xs) + (j + 1 - Ys) * (Xe - Xs) + (k + 1 - Zs) * (Xe - Xs) * (Ye - Ys);
155:           cell[7] = (i - Xs) + (j + 1 - Ys) * (Xe - Xs) + (k + 1 - Zs) * (Xe - Xs) * (Ye - Ys);
156:           if (da->elementtype == DMDA_ELEMENT_P1) {
157:             for (c = 0; c < ns * nn; c++) da->e[cnt++] = cell[split[c]];
158:           }
159:           if (da->elementtype == DMDA_ELEMENT_Q1) {
160:             for (c = 0; c < ns * nn; c++) da->e[cnt++] = cell[c];
161:           }
162:         }
163:       }
164:     }

166:     corners[0] = (xs - Xs) + (ys - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys);
167:     corners[1] = (xe - 1 - Xs) + (ys - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys);
168:     corners[2] = (xs - Xs) + (ye - 1 - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys);
169:     corners[3] = (xe - 1 - Xs) + (ye - 1 - Ys) * (Xe - Xs) + (zs - Zs) * (Xe - Xs) * (Ye - Ys);
170:     corners[4] = (xs - Xs) + (ys - Ys) * (Xe - Xs) + (ze - 1 - Zs) * (Xe - Xs) * (Ye - Ys);
171:     corners[5] = (xe - 1 - Xs) + (ys - Ys) * (Xe - Xs) + (ze - 1 - Zs) * (Xe - Xs) * (Ye - Ys);
172:     corners[6] = (xs - Xs) + (ye - 1 - Ys) * (Xe - Xs) + (ze - 1 - Zs) * (Xe - Xs) * (Ye - Ys);
173:     corners[7] = (xe - 1 - Xs) + (ye - 1 - Ys) * (Xe - Xs) + (ze - 1 - Zs) * (Xe - Xs) * (Ye - Ys);
174:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 8, corners, PETSC_COPY_VALUES, &da->ecorners));
175:   }
176:   *nel = da->ne;
177:   *nen = da->nen;
178:   *e   = da->e;
179:   PetscFunctionReturn(PETSC_SUCCESS);
180: }

182: /*@
183:    DMDAGetElementsCorners - Returns the global (x,y,z) indices of the lower left
184:    corner of the non-overlapping decomposition identified by `DMDAGetElements()`

186:     Not Collective

188:    Input Parameter:
189: .     da - the `DMDA` object

191:    Output Parameters:
192: +     gx - the x index
193: .     gy - the y index
194: -     gz - the z index

196:    Level: intermediate

198: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`
199: @*/
200: PetscErrorCode DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
201: {
202:   PetscInt  xs, Xs;
203:   PetscInt  ys, Ys;
204:   PetscInt  zs, Zs;
205:   PetscBool isda;

207:   PetscFunctionBegin;
212:   PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda));
213:   PetscCheck(isda, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)da)->type_name);
214:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, NULL, NULL, NULL));
215:   PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
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 `DMDA` 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:    Note:
241:     It returns the same number of elements, irrespective of the `DMDAElementType`

243: .seealso: `DM`, `DMDA`, `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;

253:   PetscFunctionBegin;
258:   PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda));
259:   PetscCheck(isda, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)da)->type_name);
260:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xe, &ye, &ze));
261:   PetscCall(DMDAGetGhostCorners(da, &Xs, &Ys, &Zs, NULL, NULL, NULL));
262:   xe += xs;
263:   if (xs != Xs) xs -= 1;
264:   ye += ys;
265:   if (ys != Ys) ys -= 1;
266:   ze += zs;
267:   if (zs != Zs) zs -= 1;
268:   if (mx) *mx = 0;
269:   if (my) *my = 0;
270:   if (mz) *mz = 0;
271:   PetscCall(DMGetDimension(da, &dim));
272:   switch (dim) {
273:   case 3:
274:     if (mz) *mz = ze - zs - 1; /* fall through */
275:   case 2:
276:     if (my) *my = ye - ys - 1; /* fall through */
277:   case 1:
278:     if (mx) *mx = xe - xs - 1;
279:     break;
280:   }
281:   PetscFunctionReturn(PETSC_SUCCESS);
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 Parameter:
293: .     etype - the element type, currently either `DMDA_ELEMENT_P1` or `DMDA_ELEMENT_Q1`

295:    Level: intermediate

297: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDAGetElementType()`, `DMDAGetElements()`, `DMDARestoreElements()`
298: @*/
299: PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype)
300: {
301:   DM_DA    *dd = (DM_DA *)da->data;
302:   PetscBool isda;

304:   PetscFunctionBegin;
307:   PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda));
308:   if (!isda) PetscFunctionReturn(PETSC_SUCCESS);
309:   if (dd->elementtype != etype) {
310:     PetscCall(PetscFree(dd->e));
311:     PetscCall(ISDestroy(&dd->ecorners));

313:     dd->elementtype = etype;
314:     dd->ne          = 0;
315:     dd->nen         = 0;
316:     dd->e           = NULL;
317:   }
318:   PetscFunctionReturn(PETSC_SUCCESS);
319: }

321: /*@
322:       DMDAGetElementType - Gets the element type to be returned by `DMDAGetElements()`

324:     Not Collective

326:    Input Parameter:
327: .     da - the `DMDA` object

329:    Output Parameter:
330: .     etype - the element type, currently either `DMDA_ELEMENT_P1` or `DMDA_ELEMENT_Q1`

332:    Level: intermediate

334: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`, `DMDARestoreElements()`
335: @*/
336: PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype)
337: {
338:   DM_DA    *dd = (DM_DA *)da->data;
339:   PetscBool isda;

341:   PetscFunctionBegin;
344:   PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda));
345:   PetscCheck(isda, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)da)->type_name);
346:   *etype = dd->elementtype;
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: /*@C
351:       DMDAGetElements - Gets an array containing the indices (in local coordinates)
352:                  of all the local elements

354:     Not Collective; No Fortran Support

356:    Input Parameter:
357: .     dm - the `DMDA` 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: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `VecSetValuesLocal()`, `MatSetValuesLocal()`, `DMGlobalToLocalBegin()`, `DMLocalToGlobalBegin()`
374: @*/
375: PetscErrorCode DMDAGetElements(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[])
376: {
377:   PetscInt  dim;
378:   DM_DA    *dd = (DM_DA *)dm->data;
379:   PetscBool isda;

381:   PetscFunctionBegin;
386:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
387:   PetscCheck(isda, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)dm)->type_name);
388:   PetscCheck(dd->stencil_type != DMDA_STENCIL_STAR, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DMDAGetElements() requires you use a stencil type of DMDA_STENCIL_BOX");
389:   PetscCall(DMGetDimension(dm, &dim));
390:   if (dd->e) {
391:     *nel = dd->ne;
392:     *nen = dd->nen;
393:     *e   = dd->e;
394:     PetscFunctionReturn(PETSC_SUCCESS);
395:   }
396:   if (dim == -1) {
397:     *nel = 0;
398:     *nen = 0;
399:     *e   = NULL;
400:   } else if (dim == 1) {
401:     PetscCall(DMDAGetElements_1D(dm, nel, nen, e));
402:   } else if (dim == 2) {
403:     PetscCall(DMDAGetElements_2D(dm, nel, nen, e));
404:   } else if (dim == 3) {
405:     PetscCall(DMDAGetElements_3D(dm, nel, nen, e));
406:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "DMDA dimension not 1, 2, or 3, it is %" PetscInt_FMT, dim);
407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }

410: /*@
411:       DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates)
412:                                  of the non-overlapping decomposition identified by `DMDAGetElements()`

414:     Not Collective

416:    Input Parameter:
417: .     dm - the `DMDA` object

419:    Output Parameter:
420: .     is - the index set

422:    Level: intermediate

424:    Note:
425:     Call `DMDARestoreSubdomainCornersIS()` once you have finished accessing the index set.

427: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`, `DMDARestoreElementsCornersIS()`
428: @*/
429: PetscErrorCode DMDAGetSubdomainCornersIS(DM dm, IS *is)
430: {
431:   DM_DA    *dd = (DM_DA *)dm->data;
432:   PetscBool isda;

434:   PetscFunctionBegin;
437:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
438:   PetscCheck(isda, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Not for DM type %s", ((PetscObject)dm)->type_name);
439:   PetscCheck(dd->stencil_type != DMDA_STENCIL_STAR, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DMDAGetElement() requires you use a stencil type of DMDA_STENCIL_BOX");
440:   if (!dd->ecorners) { /* compute elements if not yet done */
441:     const PetscInt *e;
442:     PetscInt        nel, nen;

444:     PetscCall(DMDAGetElements(dm, &nel, &nen, &e));
445:     PetscCall(DMDARestoreElements(dm, &nel, &nen, &e));
446:   }
447:   *is = dd->ecorners;
448:   PetscFunctionReturn(PETSC_SUCCESS);
449: }

451: /*@C
452:       DMDARestoreElements - Restores the array obtained with `DMDAGetElements()`

454:     Not Collective; No Fortran Support

456:    Input Parameters:
457: +     dm - the `DM` object
458: .     nel - number of local elements
459: .     nen - number of element nodes
460: -     e - the local indices of the elements' vertices

462:    Level: intermediate

464:    Note:
465:    This restore signals the `DMDA` object that you no longer need access to the array information.

467: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetElements()`
468: @*/
469: PetscErrorCode DMDARestoreElements(DM dm, PetscInt *nel, PetscInt *nen, const PetscInt *e[])
470: {
471:   PetscFunctionBegin;
476:   *nel = 0;
477:   *nen = -1;
478:   *e   = NULL;
479:   PetscFunctionReturn(PETSC_SUCCESS);
480: }

482: /*@
483:       DMDARestoreSubdomainCornersIS - Restores the `IS` obtained with `DMDAGetSubdomainCornersIS()`

485:     Not Collective

487:    Input Parameters:
488: +     dm - the `DM` object
489: -     is - the index set

491:    Level: intermediate

493: .seealso: `DM`, `DMDA`, `DMDAElementType`, `DMDASetElementType()`, `DMDAGetSubdomainCornersIS()`
494: @*/
495: PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm, IS *is)
496: {
497:   PetscFunctionBegin;
500:   *is = NULL;
501:   PetscFunctionReturn(PETSC_SUCCESS);
502: }