Actual source code: dagetelem.c

petsc-3.8.3 2017-12-09
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:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
 14:     DMDAGetCorners(dm,&xs,0,0,&xe,0,0);
 15:     DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);
 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:   }
 24:   *nel = da->ne;
 25:   *nen = 2;
 26:   *e   = da->e;
 27:   return(0);
 28: }

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

 41:   if (da->elementtype == DMDA_ELEMENT_P1) {nn=3;}
 42:   if (da->elementtype == DMDA_ELEMENT_Q1) {nn=4;}
 43:   if (!da->e) {
 44:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
 45:     if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;}
 46:     if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
 47:     DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);
 48:     DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);
 49:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
 50:     ye    += ys; Ye += Ys; if (ys != Ys) ys -= 1;
 51:     da->ne = ns*(xe - xs - 1)*(ye - ys - 1);
 52:     PetscMalloc1(1 + nn*da->ne,&da->e);
 53:     for (j=ys; j<ye-1; j++) {
 54:       for (i=xs; i<xe-1; i++) {
 55:         cell[0] = (i-Xs  ) + (j-Ys  )*(Xe-Xs);
 56:         cell[1] = (i-Xs+1) + (j-Ys  )*(Xe-Xs);
 57:         cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs);
 58:         cell[3] = (i-Xs  ) + (j-Ys+1)*(Xe-Xs);
 59:         if (da->elementtype == DMDA_ELEMENT_P1) {
 60:           for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
 61:         }
 62:         if (da->elementtype == DMDA_ELEMENT_Q1) {
 63:           for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
 64:         }
 65:       }
 66:     }
 67:   }
 68:   *nel = da->ne;
 69:   *nen = nn;
 70:   *e   = da->e;
 71:   return(0);
 72: }

 74: static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
 75: {
 77:   DM_DA          *da = (DM_DA*)dm->data;
 78:   PetscInt       i,xs,xe,Xs,Xe;
 79:   PetscInt       j,ys,ye,Ys,Ye;
 80:   PetscInt       k,zs,ze,Zs,Ze;
 81:   PetscInt       cnt=0, cell[8], ns=6, nn=4;
 82:   PetscInt       c, split[] = {0,1,3,7,
 83:                                0,1,7,4,
 84:                                1,2,3,7,
 85:                                1,2,7,6,
 86:                                1,4,5,7,
 87:                                1,5,6,7};

 90:   if (da->elementtype == DMDA_ELEMENT_P1) {nn=4;}
 91:   if (da->elementtype == DMDA_ELEMENT_Q1) {nn=8;}
 92:   if (!da->e) {
 93:     if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width");
 94:     if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;}
 95:     if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;}
 96:     DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);
 97:     DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);
 98:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
 99:     ye    += ys; Ye += Ys; if (ys != Ys) ys -= 1;
100:     ze    += zs; Ze += Zs; if (zs != Zs) zs -= 1;
101:     da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1);
102:     PetscMalloc1(1 + nn*da->ne,&da->e);
103:     for (k=zs; k<ze-1; k++) {
104:       for (j=ys; j<ye-1; j++) {
105:         for (i=xs; i<xe-1; i++) {
106:           cell[0] = (i-Xs  ) + (j-Ys  )*(Xe-Xs) + (k-Zs  )*(Xe-Xs)*(Ye-Ys);
107:           cell[1] = (i-Xs+1) + (j-Ys  )*(Xe-Xs) + (k-Zs  )*(Xe-Xs)*(Ye-Ys);
108:           cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs  )*(Xe-Xs)*(Ye-Ys);
109:           cell[3] = (i-Xs  ) + (j-Ys+1)*(Xe-Xs) + (k-Zs  )*(Xe-Xs)*(Ye-Ys);
110:           cell[4] = (i-Xs  ) + (j-Ys  )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
111:           cell[5] = (i-Xs+1) + (j-Ys  )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
112:           cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
113:           cell[7] = (i-Xs  ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys);
114:           if (da->elementtype == DMDA_ELEMENT_P1) {
115:             for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]];
116:           }
117:           if (da->elementtype == DMDA_ELEMENT_Q1) {
118:             for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c];
119:           }
120:         }
121:       }
122:     }
123:   }
124:   *nel = da->ne;
125:   *nen = nn;
126:   *e   = da->e;
127:   return(0);
128: }

130: /*@C
131:       DMDASetElementType - Sets the element type to be returned by DMDAGetElements()

133:     Not Collective

135:    Input Parameter:
136: .     da - the DMDA object

138:    Output Parameters:
139: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

141:    Level: intermediate

143: .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements()
144: @*/
145: PetscErrorCode  DMDASetElementType(DM da, DMDAElementType etype)
146: {
147:   DM_DA          *dd = (DM_DA*)da->data;

153:   if (dd->elementtype != etype) {
154:     PetscFree(dd->e);

156:     dd->elementtype = etype;
157:     dd->ne          = 0;
158:     dd->e           = NULL;
159:   }
160:   return(0);
161: }

163: /*@C
164:       DMDAGetElementType - Gets the element type to be returned by DMDAGetElements()

166:     Not Collective

168:    Input Parameter:
169: .     da - the DMDA object

171:    Output Parameters:
172: .     etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1

174:    Level: intermediate

176: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
177: @*/
178: PetscErrorCode  DMDAGetElementType(DM da, DMDAElementType *etype)
179: {
180:   DM_DA *dd = (DM_DA*)da->data;

185:   *etype = dd->elementtype;
186:   return(0);
187: }

189: /*@C
190:       DMDAGetElements - Gets an array containing the indices (in local coordinates)
191:                  of all the local elements

193:     Not Collective

195:    Input Parameter:
196: .     dm - the DM object

198:    Output Parameters:
199: +     nel - number of local elements
200: .     nen - number of element nodes
201: -     e - the local indices of the elements' vertices

203:    Level: intermediate

205:    Notes:
206:      Call DMDARestoreElements() once you have finished accessing the elements.

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

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

212: .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
213: @*/
214: PetscErrorCode  DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
215: {
216:   PetscInt       dim;
218:   DM_DA          *dd = (DM_DA*)dm->data;

221:   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");
222:   DMGetDimension(dm, &dim);
223:   if (dim==-1) {
224:     *nel = 0; *nen = 0; *e = NULL;
225:   } else if (dim==1) {
226:     DMDAGetElements_1D(dm,nel,nen,e);
227:   } else if (dim==2) {
228:     DMDAGetElements_2D(dm,nel,nen,e);
229:   } else if (dim==3) {
230:     DMDAGetElements_3D(dm,nel,nen,e);
231:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim);
232:   return(0);
233: }

235: /*@C
236:       DMDARestoreElements - Restores the array obtained with DMDAGetElements()

238:     Not Collective

240:    Input Parameter:
241: +     dm - the DM object
242: .     nel - number of local elements
243: .     nen - number of element nodes
244: -     e - the local indices of the elements' vertices

246:    Level: intermediate

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

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

252: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
253: @*/
254: PetscErrorCode  DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
255: {
261:   *nel = 0;
262:   *nen = -1;
263:   *e = NULL;
264:   return(0);
265: }