Actual source code: dagetelem.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  2: #include <petsc-private/dmdaimpl.h>     /*I  "petscdmda.h"   I*/

  6: static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
  7: {
  9:   DM_DA          *da = (DM_DA*)dm->data;
 10:   PetscInt       i,xs,xe,Xs,Xe;
 11:   PetscInt       cnt=0;

 14:   if (!da->e) {
 15:     DMDAGetCorners(dm,&xs,0,0,&xe,0,0);
 16:     DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);
 17:     xe    += xs; Xe += Xs; if (xs != Xs) xs -= 1;
 18:     da->ne = 1*(xe - xs - 1);
 19:     PetscMalloc1((1 + 2*da->ne),&da->e);
 20:     for (i=xs; i<xe-1; i++) {
 21:       da->e[cnt++] = (i-Xs);
 22:       da->e[cnt++] = (i-Xs+1);
 23:     }
 24:   }
 25:   *nel = da->ne;
 26:   *nen = 2;
 27:   *e   = da->e;
 28:   return(0);
 29: }

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

 44:   if (!da->e) {
 45:     if (da->elementtype == DMDA_ELEMENT_P1) {ns=2; nn=3;}
 46:     if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1; nn=4;}
 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: }

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

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

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

132:     Not Collective

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

137:    Output Parameters:
138: .     etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1

140:    Level: intermediate

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

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

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

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

167:     Not Collective

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

172:    Output Parameters:
173: .     etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1

175:    Level: intermediate

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

188:   *etype = dd->elementtype;
189:   return(0);
190: }

192: /*@C
193:       DMDAGetElements - Gets an array containing the indices (in local coordinates)
194:                  of all the local elements

196:     Not Collective

198:    Input Parameter:
199: .     dm - the DM object

201:    Output Parameters:
202: +     nel - number of local elements
203: .     nen - number of element nodes
204: -     e - the local indices of the elements' vertices

206:    Level: intermediate

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

218:   if (da->dim==-1) {
219:     *nel = 0; *nen = 0; *e = NULL;
220:   } else if (da->dim==1) {
221:     DMDAGetElements_1D(dm,nel,nen,e);
222:   } else if (da->dim==2) {
223:     DMDAGetElements_2D(dm,nel,nen,e);
224:   } else if (da->dim==3) {
225:     DMDAGetElements_3D(dm,nel,nen,e);
226:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",da->dim);
227:   return(0);
228: }

232: /*@C
233:       DMDARestoreElements - Returns an array containing the indices (in local coordinates)
234:                  of all the local elements obtained with DMDAGetElements()

236:     Not Collective

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

244:    Level: intermediate

246: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
247: @*/
248: PetscErrorCode  DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
249: {
255:   /* XXX */
256:   return(0);
257: }