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; Xe += Xs; if (xs != Xs) xs -= 1;
 17:     da->ne = 1*(xe - xs - 1);
 18:     PetscMalloc1(1 + 2*da->ne,&da->e);
 19:     for (i=xs; i<xe-1; i++) {
 20:       da->e[cnt++] = (i-Xs);
 21:       da->e[cnt++] = (i-Xs+1);
 22:     }
 23:     da->nen = 2;

 25:     corners[0] = (xs  -Xs);
 26:     corners[1] = (xe-1-Xs);
 27:     ISCreateGeneral(PETSC_COMM_SELF,2,corners,PETSC_COPY_VALUES,&da->ecorners);
 28:   }
 29:   *nel = da->ne;
 30:   *nen = da->nen;
 31:   *e   = da->e;
 32:   return 0;
 33: }

 35: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
 36: {
 37:   DM_DA          *da = (DM_DA*)dm->data;
 38:   PetscInt       i,xs,xe,Xs,Xe;
 39:   PetscInt       j,ys,ye,Ys,Ye;
 40:   PetscInt       cnt=0, cell[4], ns=2;
 41:   PetscInt       c, split[] = {0,1,3,
 42:                                2,3,1};

 44:   if (!da->e) {
 45:     PetscInt corners[4],nn = 0;


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

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

 84:     corners[0] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs);
 85:     corners[1] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs);
 86:     corners[2] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs);
 87:     corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs);
 88:     ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);
 89:   }
 90:   *nel = da->ne;
 91:   *nen = da->nen;
 92:   *e   = da->e;
 93:   return 0;
 94: }

 96: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
 97: {
 98:   DM_DA          *da = (DM_DA*)dm->data;
 99:   PetscInt       i,xs,xe,Xs,Xe;
100:   PetscInt       j,ys,ye,Ys,Ye;
101:   PetscInt       k,zs,ze,Zs,Ze;
102:   PetscInt       cnt=0, cell[8], ns=6;
103:   PetscInt       c, split[] = {0,1,3,7,
104:                                0,1,7,4,
105:                                1,2,3,7,
106:                                1,2,7,6,
107:                                1,4,5,7,
108:                                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; Xe += Xs; if (xs != Xs) xs -= 1;
132:     ye    += ys; Ye += Ys; if (ys != Ys) ys -= 1;
133:     ze    += zs; Ze += Zs; if (zs != Zs) zs -= 1;
134:     da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
135:     PetscMalloc1(1 + nn*da->ne,&da->e);
136:     for (k=zs; k<ze-1; k++) {
137:       for (j=ys; j<ye-1; j++) {
138:         for (i=xs; i<xe-1; i++) {
139:           cell[0] = (i  -Xs) + (j  -Ys)*(Xe-Xs) + (k  -Zs)*(Xe-Xs)*(Ye-Ys);
140:           cell[1] = (i+1-Xs) + (j  -Ys)*(Xe-Xs) + (k  -Zs)*(Xe-Xs)*(Ye-Ys);
141:           cell[2] = (i+1-Xs) + (j+1-Ys)*(Xe-Xs) + (k  -Zs)*(Xe-Xs)*(Ye-Ys);
142:           cell[3] = (i  -Xs) + (j+1-Ys)*(Xe-Xs) + (k  -Zs)*(Xe-Xs)*(Ye-Ys);
143:           cell[4] = (i  -Xs) + (j  -Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
144:           cell[5] = (i+1-Xs) + (j  -Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
145:           cell[6] = (i+1-Xs) + (j+1-Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
146:           cell[7] = (i  -Xs) + (j+1-Ys)*(Xe-Xs) + (k+1-Zs)*(Xe-Xs)*(Ye-Ys);
147:           if (da->elementtype == DMDA_ELEMENT_P1) {
148:             for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
149:           }
150:           if (da->elementtype == DMDA_ELEMENT_Q1) {
151:             for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
152:           }
153:         }
154:       }
155:     }

157:     corners[0] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs) + (zs-  Zs)*(Xe-Xs)*(Ye-Ys);
158:     corners[1] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs) + (zs-  Zs)*(Xe-Xs)*(Ye-Ys);
159:     corners[2] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-  Zs)*(Xe-Xs)*(Ye-Ys);
160:     corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-  Zs)*(Xe-Xs)*(Ye-Ys);
161:     corners[4] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
162:     corners[5] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
163:     corners[6] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
164:     corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
165:     ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);
166:   }
167:   *nel = da->ne;
168:   *nen = da->nen;
169:   *e   = da->e;
170:   return 0;
171: }

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

177:     Not Collective

179:    Input Parameter:
180: .     da - the DM object

182:    Output Parameters:
183: +     gx - the x index
184: .     gy - the y index
185: -     gz - the z index

187:    Level: intermediate

189:    Notes:

191: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
192: @*/
193: PetscErrorCode  DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
194: {
195:   PetscInt       xs,Xs;
196:   PetscInt       ys,Ys;
197:   PetscInt       zs,Zs;
198:   PetscBool      isda;

204:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
206:   DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);
207:   DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
208:   if (xs != Xs) xs -= 1;
209:   if (ys != Ys) ys -= 1;
210:   if (zs != Zs) zs -= 1;
211:   if (gx) *gx  = xs;
212:   if (gy) *gy  = ys;
213:   if (gz) *gz  = zs;
214:   return 0;
215: }

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

220:     Not Collective

222:    Input Parameter:
223: .     da - the DM object

225:    Output Parameters:
226: +     mx - number of local elements in x-direction
227: .     my - number of local elements in y-direction
228: -     mz - number of local elements in z-direction

230:    Level: intermediate

232:    Notes:
233:     It returns the same number of elements, irrespective of the DMDAElementType

235: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements
236: @*/
237: PetscErrorCode  DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
238: {
239:   PetscInt       xs,xe,Xs;
240:   PetscInt       ys,ye,Ys;
241:   PetscInt       zs,ze,Zs;
242:   PetscInt       dim;
243:   PetscBool      isda;

249:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
251:   DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);
252:   DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
253:   xe  += xs; if (xs != Xs) xs -= 1;
254:   ye  += ys; if (ys != Ys) ys -= 1;
255:   ze  += zs; if (zs != Zs) zs -= 1;
256:   if (mx) *mx  = 0;
257:   if (my) *my  = 0;
258:   if (mz) *mz  = 0;
259:   DMGetDimension(da,&dim);
260:   switch (dim) {
261:   case 3:
262:     if (mz) *mz = ze - zs - 1;
263:   case 2:
264:     if (my) *my = ye - ys - 1;
265:   case 1:
266:     if (mx) *mx = xe - xs - 1;
267:     break;
268:   }
269:   return 0;
270: }

272: /*@
273:       DMDASetElementType - Sets the element type to be returned by DMDAGetElements()

275:     Not Collective

277:    Input Parameter:
278: .     da - the DMDA object

280:    Output Parameters:
281: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

283:    Level: intermediate

285: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
286: @*/
287: PetscErrorCode  DMDASetElementType(DM da, DMDAElementType etype)
288: {
289:   DM_DA          *dd = (DM_DA*)da->data;
290:   PetscBool      isda;

294:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
295:   if (!isda) return 0;
296:   if (dd->elementtype != etype) {
297:     PetscFree(dd->e);
298:     ISDestroy(&dd->ecorners);

300:     dd->elementtype = etype;
301:     dd->ne          = 0;
302:     dd->nen         = 0;
303:     dd->e           = NULL;
304:   }
305:   return 0;
306: }

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

311:     Not Collective

313:    Input Parameter:
314: .     da - the DMDA object

316:    Output Parameters:
317: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

319:    Level: intermediate

321: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
322: @*/
323: PetscErrorCode  DMDAGetElementType(DM da, DMDAElementType *etype)
324: {
325:   DM_DA          *dd = (DM_DA*)da->data;
326:   PetscBool      isda;

330:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
332:   *etype = dd->elementtype;
333:   return 0;
334: }

336: /*@C
337:       DMDAGetElements - Gets an array containing the indices (in local coordinates)
338:                  of all the local elements

340:     Not Collective

342:    Input Parameter:
343: .     dm - the DM object

345:    Output Parameters:
346: +     nel - number of local elements
347: .     nen - number of element nodes
348: -     e - the local indices of the elements' vertices

350:    Level: intermediate

352:    Notes:
353:      Call DMDARestoreElements() once you have finished accessing the elements.

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

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

359:      Not supported in Fortran

361: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
362: @*/
363: PetscErrorCode  DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
364: {
365:   PetscInt       dim;
366:   DM_DA          *dd = (DM_DA*)dm->data;
367:   PetscBool      isda;

373:   PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
376:   DMGetDimension(dm, &dim);
377:   if (dd->e) {
378:     *nel = dd->ne;
379:     *nen = dd->nen;
380:     *e   = dd->e;
381:     return 0;
382:   }
383:   if (dim==-1) {
384:     *nel = 0; *nen = 0; *e = NULL;
385:   } else if (dim==1) {
386:     DMDAGetElements_1D(dm,nel,nen,e);
387:   } else if (dim==2) {
388:     DMDAGetElements_2D(dm,nel,nen,e);
389:   } else if (dim==3) {
390:     DMDAGetElements_3D(dm,nel,nen,e);
391:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D",dim);
392:   return 0;
393: }

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

399:     Not Collective

401:    Input Parameter:
402: .     dm - the DM object

404:    Output Parameters:
405: .     is - the index set

407:    Level: intermediate

409:    Notes:
410:     Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set.

412: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS()
413: @*/
414: PetscErrorCode  DMDAGetSubdomainCornersIS(DM dm,IS *is)
415: {
416:   DM_DA          *dd = (DM_DA*)dm->data;
417:   PetscBool      isda;

421:   PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
424:   if (!dd->ecorners) { /* compute elements if not yet done */
425:     const PetscInt *e;
426:     PetscInt       nel,nen;

428:     DMDAGetElements(dm,&nel,&nen,&e);
429:     DMDARestoreElements(dm,&nel,&nen,&e);
430:   }
431:   *is = dd->ecorners;
432:   return 0;
433: }

435: /*@C
436:       DMDARestoreElements - Restores the array obtained with DMDAGetElements()

438:     Not Collective

440:    Input Parameters:
441: +     dm - the DM object
442: .     nel - number of local elements
443: .     nen - number of element nodes
444: -     e - the local indices of the elements' vertices

446:    Level: intermediate

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

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

452:          Not supported in Fortran

454: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
455: @*/
456: PetscErrorCode  DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
457: {
462:   *nel = 0;
463:   *nen = -1;
464:   *e = NULL;
465:   return 0;
466: }

468: /*@
469:       DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS()

471:     Not Collective

473:    Input Parameters:
474: +     dm - the DM object
475: -     is - the index set

477:    Level: intermediate

479:    Note:

481: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS()
482: @*/
483: PetscErrorCode  DMDARestoreSubdomainCornersIS(DM dm,IS *is)
484: {
487:   *is = NULL;
488:   return 0;
489: }