Actual source code: dagetelem.c

petsc-3.11.1 2019-04-17
Report Typos and Errors

  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:     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:       break;
 63:     }
 64:     nn = da->nen;

 66:     if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
 67:     if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
 68:     DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);
 69:     DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);
 70:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
 71:     ye    += ys; Ye += Ys; 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: {
104:   DM_DA          *da = (DM_DA*)dm->data;
105:   PetscInt       i,xs,xe,Xs,Xe;
106:   PetscInt       j,ys,ye,Ys,Ye;
107:   PetscInt       k,zs,ze,Zs,Ze;
108:   PetscInt       cnt=0, cell[8], ns=6;
109:   PetscInt       c, split[] = {0,1,3,7,
110:                                0,1,7,4,
111:                                1,2,3,7,
112:                                1,2,7,6,
113:                                1,4,5,7,
114:                                1,5,6,7};

117:   if (!da->e) {
118:     PetscInt corners[8],nn = 0;

120:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");

122:     switch (da->elementtype) {
123:     case DMDA_ELEMENT_Q1:
124:       da->nen = 8;
125:       break;
126:     case DMDA_ELEMENT_P1:
127:       da->nen = 4;
128:       break;
129:     default:
130:       SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unknown element type %d",da->elementtype);
131:       break;
132:     }
133:     nn = da->nen;

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

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

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

185:     Not Collective

187:    Input Parameter:
188: .     da - the DM object

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

195:    Level: intermediate

197:    Notes:

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

214:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
215:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
216:   DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);
217:   DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
218:   if (xs != Xs) xs -= 1;
219:   if (ys != Ys) ys -= 1;
220:   if (zs != Zs) zs -= 1;
221:   if (gx) *gx  = xs;
222:   if (gy) *gy  = ys;
223:   if (gz) *gz  = zs;
224:   return(0);
225: }

227: /*@
228:       DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by DMDAGetElements()

230:     Not Collective

232:    Input Parameter:
233: .     da - the DM object

235:    Output Parameters:
236: +     mx - number of local elements in x-direction
237: .     my - number of local elements in y-direction
238: -     mz - number of local elements in z-direction

240:    Level: intermediate

242:    Notes:
243:     It returns the same number of elements, irrespective of the DMDAElementType

245: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements
246: @*/
247: PetscErrorCode  DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
248: {
249:   PetscInt       xs,xe,Xs;
250:   PetscInt       ys,ye,Ys;
251:   PetscInt       zs,ze,Zs;
252:   PetscInt       dim;
253:   PetscBool      isda;

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

295:    Level: intermediate

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

308:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
309:   if (!isda) return(0);
310:   if (dd->elementtype != etype) {
311:     PetscFree(dd->e);
312:     ISDestroy(&dd->ecorners);

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

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

325:     Not Collective

327:    Input Parameter:
328: .     da - the DMDA object

330:    Output Parameters:
331: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

333:    Level: intermediate

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

346:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
347:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
348:   *etype = dd->elementtype;
349:   return(0);
350: }

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

356:     Not Collective

358:    Input Parameter:
359: .     dm - the DM object

361:    Output Parameters:
362: +     nel - number of local elements
363: .     nen - number of element nodes
364: -     e - the local indices of the elements' vertices

366:    Level: intermediate

368:    Notes:
369:      Call DMDARestoreElements() once you have finished accessing the elements.

371:      Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes.

373:      If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result.

375:      Not supported in Fortran

377: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
378: @*/
379: PetscErrorCode  DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
380: {
381:   PetscInt       dim;
383:   DM_DA          *dd = (DM_DA*)dm->data;
384:   PetscBool      isda;

391:   PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
392:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
393:   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");
394:   DMGetDimension(dm, &dim);
395:   if (dd->e) {
396:     *nel = dd->ne;
397:     *nen = dd->nen;
398:     *e   = dd->e;
399:     return(0);
400:   }
401:   if (dim==-1) {
402:     *nel = 0; *nen = 0; *e = NULL;
403:   } else if (dim==1) {
404:     DMDAGetElements_1D(dm,nel,nen,e);
405:   } else if (dim==2) {
406:     DMDAGetElements_2D(dm,nel,nen,e);
407:   } else if (dim==3) {
408:     DMDAGetElements_3D(dm,nel,nen,e);
409:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
410:   return(0);
411: }

413: /*@
414:       DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates)
415:                                  of the non-overlapping decomposition identified by DMDAGetElements

417:     Not Collective

419:    Input Parameter:
420: .     dm - the DM object

422:    Output Parameters:
423: .     is - the index set

425:    Level: intermediate

427:    Notes:
428:     Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set.

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

441:   PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
442:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
443:   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");
444:   if (!dd->ecorners) { /* compute elements if not yet done */
445:     const PetscInt *e;
446:     PetscInt       nel,nen;

448:     DMDAGetElements(dm,&nel,&nen,&e);
449:     DMDARestoreElements(dm,&nel,&nen,&e);
450:   }
451:   *is = dd->ecorners;
452:   return(0);
453: }

455: /*@C
456:       DMDARestoreElements - Restores the array obtained with DMDAGetElements()

458:     Not Collective

460:    Input Parameter:
461: +     dm - the DM object
462: .     nel - number of local elements
463: .     nen - number of element nodes
464: -     e - the local indices of the elements' vertices

466:    Level: intermediate

468:    Note: You should not access these values after you have called this routine.

470:          This restore signals the DMDA object that you no longer need access to the array information.

472:          Not supported in Fortran

474: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
475: @*/
476: PetscErrorCode  DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
477: {
483:   *nel = 0;
484:   *nen = -1;
485:   *e = NULL;
486:   return(0);
487: }

489: /*@
490:       DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS()

492:     Not Collective

494:    Input Parameter:
495: +     dm - the DM object
496: -     is - the index set

498:    Level: intermediate

500:    Note:

502: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS()
503: @*/
504: PetscErrorCode  DMDARestoreSubdomainCornersIS(DM dm,IS *is)
505: {
509:   *is = NULL;
510:   return(0);
511: }