Actual source code: dagetelem.c

petsc-3.9.0 2018-04-07
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:     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 = 2;
 31:   *e   = da->e;
 32:   return(0);
 33: }

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

 46:   if (da->elementtype == DMDA_ELEMENT_P1) {nn=3;}
 47:   if (da->elementtype == DMDA_ELEMENT_Q1) {nn=4;}
 48:   if (!da->e) {
 49:     PetscInt corners[4];

 51:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
 52:     if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
 53:     if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
 54:     DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);
 55:     DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);
 56:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
 57:     ye    += ys; Ye += Ys; if (ys != Ys) ys -= 1;
 58:     da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
 59:     PetscMalloc1(1 + nn*da->ne,&da->e);
 60:     for (j=ys; j<ye-1; j++) {
 61:       for (i=xs; i<xe-1; i++) {
 62:         cell[0] = (i-Xs  ) + (j-Ys  )*(Xe-Xs);
 63:         cell[1] = (i-Xs+1) + (j-Ys  )*(Xe-Xs);
 64:         cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
 65:         cell[3] = (i-Xs  ) + (j-Ys+1)*(Xe-Xs);
 66:         if (da->elementtype == DMDA_ELEMENT_P1) {
 67:           for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
 68:         }
 69:         if (da->elementtype == DMDA_ELEMENT_Q1) {
 70:           for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
 71:         }
 72:       }
 73:     }
 74:     corners[0] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs);
 75:     corners[1] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs);
 76:     corners[2] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs);
 77:     corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs);
 78:     ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);
 79:   }
 80:   *nel = da->ne;
 81:   *nen = nn;
 82:   *e   = da->e;
 83:   return(0);
 84: }

 86: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
 87: {
 89:   DM_DA          *da = (DM_DA*)dm->data;
 90:   PetscInt       i,xs,xe,Xs,Xe;
 91:   PetscInt       j,ys,ye,Ys,Ye;
 92:   PetscInt       k,zs,ze,Zs,Ze;
 93:   PetscInt       cnt=0, cell[8], ns=6, nn=4;
 94:   PetscInt       c, split[] = {0,1,3,7,
 95:                                0,1,7,4,
 96:                                1,2,3,7,
 97:                                1,2,7,6,
 98:                                1,4,5,7,
 99:                                1,5,6,7};

102:   if (da->elementtype == DMDA_ELEMENT_P1) {nn=4;}
103:   if (da->elementtype == DMDA_ELEMENT_Q1) {nn=8;}
104:   if (!da->e) {
105:     PetscInt corners[8];

107:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
108:     if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;}
109:     if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
110:     DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
111:     DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
112:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
113:     ye    += ys; Ye += Ys; if (ys != Ys) ys -= 1;
114:     ze    += zs; Ze += Zs; if (zs != Zs) zs -= 1;
115:     da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
116:     PetscMalloc1(1 + nn*da->ne,&da->e);
117:     for (k=zs; k<ze-1; k++) {
118:       for (j=ys; j<ye-1; j++) {
119:         for (i=xs; i<xe-1; i++) {
120:           cell[0] = (i-Xs  ) + (j-Ys  )*(Xe-Xs) + (k-Zs  )*(Xe-Xs)*(Ye-Ys);
121:           cell[1] = (i-Xs+1) + (j-Ys  )*(Xe-Xs) + (k-Zs  )*(Xe-Xs)*(Ye-Ys);
122:           cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs  )*(Xe-Xs)*(Ye-Ys);
123:           cell[3] = (i-Xs  ) + (j-Ys+1)*(Xe-Xs) + (k-Zs  )*(Xe-Xs)*(Ye-Ys);
124:           cell[4] = (i-Xs  ) + (j-Ys  )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
125:           cell[5] = (i-Xs+1) + (j-Ys  )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
126:           cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
127:           cell[7] = (i-Xs  ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
128:           if (da->elementtype == DMDA_ELEMENT_P1) {
129:             for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
130:           }
131:           if (da->elementtype == DMDA_ELEMENT_Q1) {
132:             for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
133:           }
134:         }
135:       }
136:     }
137:     corners[0] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs) + (zs-Zs  )*(Xe-Xs)*(Ye-Ys);
138:     corners[1] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs) + (zs-Zs  )*(Xe-Xs)*(Ye-Ys);
139:     corners[2] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs  )*(Xe-Xs)*(Ye-Ys);
140:     corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs  )*(Xe-Xs)*(Ye-Ys);
141:     corners[4] = (xs  -Xs) + (ys  -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
142:     corners[5] = (xe-1-Xs) + (ys  -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
143:     corners[6] = (xs  -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
144:     corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys);
145:     ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);
146:   }
147:   *nel = da->ne;
148:   *nen = nn;
149:   *e   = da->e;
150:   return(0);
151: }

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

157:     Not Collective

159:    Input Parameter:
160: .     da - the DM object

162:    Output Parameters:
163: +     gx - the x index
164: .     gy - the y index
165: -     gz - the z index

167:    Level: intermediate

169:    Notes:

171: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
172: @*/
173: PetscErrorCode  DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz)
174: {
175:   PetscInt       xs,Xs;
176:   PetscInt       ys,Ys;
177:   PetscInt       zs,Zs;
178:   PetscBool      isda;

186:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
187:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
188:   DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);
189:   DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
190:   if (xs != Xs) xs -= 1;
191:   if (ys != Ys) ys -= 1;
192:   if (zs != Zs) zs -= 1;
193:   if (gx) *gx  = xs;
194:   if (gy) *gy  = ys;
195:   if (gz) *gz  = zs;
196:   return(0);
197: }

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

202:     Not Collective

204:    Input Parameter:
205: .     da - the DM object

207:    Output Parameters:
208: +     mx - number of local elements in x-direction
209: .     my - number of local elements in y-direction
210: -     mz - number of local elements in z-direction

212:    Level: intermediate

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

216: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements
217: @*/
218: PetscErrorCode  DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz)
219: {
220:   PetscInt       xs,xe,Xs;
221:   PetscInt       ys,ye,Ys;
222:   PetscInt       zs,ze,Zs;
223:   PetscInt       dim;
224:   PetscBool      isda;

232:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
233:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
234:   DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);
235:   DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);
236:   xe  += xs; if (xs != Xs) xs -= 1;
237:   ye  += ys; if (ys != Ys) ys -= 1;
238:   ze  += zs; if (zs != Zs) zs -= 1;
239:   if (mx) *mx  = 0;
240:   if (my) *my  = 0;
241:   if (mz) *mz  = 0;
242:   DMGetDimension(da,&dim);
243:   switch (dim) {
244:   case 3:
245:     if (mz) *mz = ze - zs - 1;
246:   case 2:
247:     if (my) *my = ye - ys - 1;
248:   case 1:
249:     if (mx) *mx = xe - xs - 1;
250:     break;
251:   }
252:   return(0);
253: }

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

258:     Not Collective

260:    Input Parameter:
261: .     da - the DMDA object

263:    Output Parameters:
264: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

266:    Level: intermediate

268: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
269: @*/
270: PetscErrorCode  DMDASetElementType(DM da, DMDAElementType etype)
271: {
272:   DM_DA          *dd = (DM_DA*)da->data;
274:   PetscBool      isda;

279:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
280:   if (!isda) return(0);
281:   if (dd->elementtype != etype) {
282:     PetscFree(dd->e);
283:     ISDestroy(&dd->ecorners);

285:     dd->elementtype = etype;
286:     dd->ne          = 0;
287:     dd->e           = NULL;
288:   }
289:   return(0);
290: }

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

295:     Not Collective

297:    Input Parameter:
298: .     da - the DMDA object

300:    Output Parameters:
301: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

303:    Level: intermediate

305: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
306: @*/
307: PetscErrorCode  DMDAGetElementType(DM da, DMDAElementType *etype)
308: {
309:   DM_DA          *dd = (DM_DA*)da->data;
311:   PetscBool      isda;

316:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
317:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name);
318:   *etype = dd->elementtype;
319:   return(0);
320: }

322: /*@C
323:       DMDAGetElements - Gets an array containing the indices (in local coordinates)
324:                  of all the local elements

326:     Not Collective

328:    Input Parameter:
329: .     dm - the DM object

331:    Output Parameters:
332: +     nel - number of local elements
333: .     nen - number of element nodes
334: -     e - the local indices of the elements' vertices

336:    Level: intermediate

338:    Notes:
339:      Call DMDARestoreElements() once you have finished accessing the elements.

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

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

345:      Not supported in Fortran

347: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
348: @*/
349: PetscErrorCode  DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
350: {
351:   PetscInt       dim;
353:   DM_DA          *dd = (DM_DA*)dm->data;
354:   PetscBool      isda;

361:   PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
362:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
363:   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");
364:   DMGetDimension(dm, &dim);
365:   if (dim==-1) {
366:     *nel = 0; *nen = 0; *e = NULL;
367:   } else if (dim==1) {
368:     DMDAGetElements_1D(dm,nel,nen,e);
369:   } else if (dim==2) {
370:     DMDAGetElements_2D(dm,nel,nen,e);
371:   } else if (dim==3) {
372:     DMDAGetElements_3D(dm,nel,nen,e);
373:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
374:   return(0);
375: }

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

381:     Not Collective

383:    Input Parameter:
384: .     dm - the DM object

386:    Output Parameters:
387: .     is - the index set

389:    Level: intermediate

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

393: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS()
394: @*/
395: PetscErrorCode  DMDAGetSubdomainCornersIS(DM dm,IS *is)
396: {
398:   DM_DA          *dd = (DM_DA*)dm->data;
399:   PetscBool      isda;

404:   PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);
405:   if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name);
406:   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");
407:   if (!dd->ecorners) { /* compute elements if not yet done */
408:     const PetscInt *e;
409:     PetscInt       nel,nen;

411:     DMDAGetElements(dm,&nel,&nen,&e);
412:     DMDARestoreElements(dm,&nel,&nen,&e);
413:   }
414:   *is = dd->ecorners;
415:   return(0);
416: }

418: /*@C
419:       DMDARestoreElements - Restores the array obtained with DMDAGetElements()

421:     Not Collective

423:    Input Parameter:
424: +     dm - the DM object
425: .     nel - number of local elements
426: .     nen - number of element nodes
427: -     e - the local indices of the elements' vertices

429:    Level: intermediate

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

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

435:          Not supported in Fortran

437: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
438: @*/
439: PetscErrorCode  DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
440: {
446:   *nel = 0;
447:   *nen = -1;
448:   *e = NULL;
449:   return(0);
450: }

452: /*@
453:       DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS()

455:     Not Collective

457:    Input Parameter:
458: +     dm - the DM object
459: -     is - the index set

461:    Level: intermediate

463:    Note:

465: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS()
466: @*/
467: PetscErrorCode  DMDARestoreSubdomainCornersIS(DM dm,IS *is)
468: {
472:   *is = NULL;
473:   return(0);
474: }