Actual source code: dagetelem.c

petsc-3.6.1 2015-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:     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:   }
 26:   *nel = da->ne;
 27:   *nen = 2;
 28:   *e   = da->e;
 29:   return(0);
 30: }

 34: static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
 35: {
 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, nn=3;
 41:   PetscInt       c, split[] = {0,1,3,
 42:                                2,3,1};

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

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

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

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

135:     Not Collective

137:    Input Parameter:
138: .     da - the DMDA object

140:    Output Parameters:
141: .     etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1

143:    Level: intermediate

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

157:   if (dd->elementtype != etype) {
158:     PetscFree(dd->e);

160:     dd->elementtype = etype;
161:     dd->ne          = 0;
162:     dd->e           = NULL;
163:   }
164:   return(0);
165: }

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

170:     Not Collective

172:    Input Parameter:
173: .     da - the DMDA object

175:    Output Parameters:
176: .     etype - the element type, currently either DMDA_ELEMENT_P1 or ELEMENT_Q1

178:    Level: intermediate

180: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements()
181: @*/
184: PetscErrorCode  DMDAGetElementType(DM da, DMDAElementType *etype)
185: {
186:   DM_DA *dd = (DM_DA*)da->data;

191:   *etype = dd->elementtype;
192:   return(0);
193: }

195: /*@C
196:       DMDAGetElements - Gets an array containing the indices (in local coordinates)
197:                  of all the local elements

199:     Not Collective

201:    Input Parameter:
202: .     dm - the DM object

204:    Output Parameters:
205: +     nel - number of local elements
206: .     nen - number of element nodes
207: -     e - the local indices of the elements' vertices

209:    Level: intermediate

211: .seealso: DMDAElementType, DMDASetElementType(), DMDARestoreElements(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin()
212: @*/
215: PetscErrorCode  DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
216: {
217:   PetscInt       dim;

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

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

240:     Not Collective

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

248:    Level: intermediate

250: .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements()
251: @*/
252: PetscErrorCode  DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[])
253: {
259:   /* XXX */
260:   return(0);
261: }