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: }