Actual source code: stagstencil.c

  1: /* Functions concerning getting and setting Vec and Mat values with DMStagStencil */
  2: #include <petsc/private/dmstagimpl.h>

  4: /* Strings corresponding the the types defined in $PETSC_DIR/include/petscdmstag.h */
  5: const char *const DMStagStencilTypes[] = {"NONE","STAR","BOX","DMStagStencilType","DM_STAG_STENCIL_",NULL};

  7: /* Strings corresponding the positions in $PETSC_DIR/include/petscdmstag.h */
  8: const char * const DMStagStencilLocations[] = {"NONE","BACK_DOWN_LEFT","BACK_DOWN","BACK_DOWN_RIGHT","BACK_LEFT","BACK","BACK_RIGHT","BACK_UP_LEFT","BACK_UP","BACK_UP_RIGHT","DOWN_LEFT","DOWN","DOWN_RIGHT","LEFT","ELEMENT","RIGHT","UP_LEFT","UP","UP_RIGHT","FRONT_DOWN_LEFT","FRONT_DOWN","FRONT_DOWN_RIGHT","FRONT_LEFT","FRONT","FRONT_RIGHT","FRONT_UP_LEFT","FRONT_UP","FRONT_UP_RIGHT","DMStagStencilLocation","",NULL};
  9: /*@C
 10:   DMStagGetLocationDOF - Get number of DOF associated with a given point in a DMStag grid

 12:   Not Collective

 14:   Input Parameters:
 15: + dm - the DMStag object
 16: - loc - grid point (see DMStagStencilLocation)

 18:   Output Parameter:
 19: . dof - the number of dof (components) living at loc in dm

 21:   Level: intermediate

 23: .seealso: DMSTAG, DMStagStencilLocation, DMStagStencil, DMDAGetDof()
 24: @*/
 25: PetscErrorCode DMStagGetLocationDOF(DM dm,DMStagStencilLocation loc,PetscInt *dof)
 26: {
 27:   PetscErrorCode        ierr;
 28:   const DM_Stag * const stag = (DM_Stag*)dm->data;
 29:   PetscInt              dim;

 33:   DMGetDimension(dm,&dim);
 34:   switch (dim) {
 35:     case 1:
 36:       switch (loc) {
 37:         case DMSTAG_LEFT:
 38:         case DMSTAG_RIGHT:
 39:           *dof = stag->dof[0]; break;
 40:         case DMSTAG_ELEMENT:
 41:           *dof = stag->dof[1]; break;
 42:         default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location %s",DMStagStencilLocations[loc]);
 43:       }
 44:       break;
 45:     case 2:
 46:       switch (loc) {
 47:         case DMSTAG_DOWN_LEFT:
 48:         case DMSTAG_DOWN_RIGHT:
 49:         case DMSTAG_UP_LEFT:
 50:         case DMSTAG_UP_RIGHT:
 51:           *dof = stag->dof[0]; break;
 52:         case DMSTAG_LEFT:
 53:         case DMSTAG_RIGHT:
 54:         case DMSTAG_UP:
 55:         case DMSTAG_DOWN:
 56:           *dof = stag->dof[1]; break;
 57:         case DMSTAG_ELEMENT:
 58:           *dof = stag->dof[2]; break;
 59:         default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location %s",DMStagStencilLocations[loc]);
 60:       }
 61:       break;
 62:     case 3:
 63:       switch (loc) {
 64:         case DMSTAG_BACK_DOWN_LEFT:
 65:         case DMSTAG_BACK_DOWN_RIGHT:
 66:         case DMSTAG_BACK_UP_LEFT:
 67:         case DMSTAG_BACK_UP_RIGHT:
 68:         case DMSTAG_FRONT_DOWN_LEFT:
 69:         case DMSTAG_FRONT_DOWN_RIGHT:
 70:         case DMSTAG_FRONT_UP_LEFT:
 71:         case DMSTAG_FRONT_UP_RIGHT:
 72:           *dof = stag->dof[0]; break;
 73:         case DMSTAG_BACK_DOWN:
 74:         case DMSTAG_BACK_LEFT:
 75:         case DMSTAG_BACK_RIGHT:
 76:         case DMSTAG_BACK_UP:
 77:         case DMSTAG_DOWN_LEFT:
 78:         case DMSTAG_DOWN_RIGHT:
 79:         case DMSTAG_UP_LEFT:
 80:         case DMSTAG_UP_RIGHT:
 81:         case DMSTAG_FRONT_DOWN:
 82:         case DMSTAG_FRONT_LEFT:
 83:         case DMSTAG_FRONT_RIGHT:
 84:         case DMSTAG_FRONT_UP:
 85:           *dof = stag->dof[1]; break;
 86:         case DMSTAG_LEFT:
 87:         case DMSTAG_RIGHT:
 88:         case DMSTAG_DOWN:
 89:         case DMSTAG_UP:
 90:         case DMSTAG_BACK:
 91:         case DMSTAG_FRONT:
 92:           *dof = stag->dof[2]; break;
 93:         case DMSTAG_ELEMENT:
 94:           *dof = stag->dof[3]; break;
 95:         default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location %s",DMStagStencilLocations[loc]);
 96:       }
 97:       break;
 98:     default : SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dimension %D",dim);
 99:   }
100:   return(0);
101: }

103: /*@C
104:   DMStagMatGetValuesStencil - retrieve local matrix entries using grid indexing

106:   Not Collective

108:   Input Parameters:
109: + dm - the DMStag object
110: . mat - the matrix
111: . nRow - number of rows
112: . posRow - grid locations (including components) of rows
113: . nCol - number of columns
114: - posCol - grid locations (including components) of columns

116:   Output Parameter:
117: . val - logically two-dimensional array of values

119:   Level: advanced

121: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
122: @*/
123: PetscErrorCode DMStagMatGetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,PetscScalar *val)
124: {
126:   PetscInt       dim;
127:   PetscInt       *ir,*ic;

132:   DMGetDimension(dm,&dim);
133:   PetscMalloc2(nRow,&ir,nCol,&ic);
134:   DMStagStencilToIndexLocal(dm,dim,nRow,posRow,ir);
135:   DMStagStencilToIndexLocal(dm,dim,nCol,posCol,ic);
136:   MatGetValuesLocal(mat,nRow,ir,nCol,ic,val);
137:   PetscFree2(ir,ic);
138:   return(0);
139: }

141: /*@C
142:   DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing

144:   Not Collective

146:   Input Parameters:
147: + dm - the DMStag object
148: . mat - the matrix
149: . nRow - number of rows
150: . posRow - grid locations (including components) of rows
151: . nCol - number of columns
152: . posCol - grid locations (including components) of columns
153: . val - logically two-dimensional array of values
154: - insertmode - INSERT_VALUES or ADD_VALUES

156:   Notes:
157:   See notes for MatSetValuesStencil()

159:   Level: intermediate

161: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatGetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
162: @*/
163: PetscErrorCode DMStagMatSetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,const PetscScalar *val,InsertMode insertMode)
164: {
166:   PetscInt       dim;
167:   PetscInt       *ir,*ic;

172:   DMGetDimension(dm,&dim);
173:   PetscMalloc2(nRow,&ir,nCol,&ic);
174:   DMStagStencilToIndexLocal(dm,dim,nRow,posRow,ir);
175:   DMStagStencilToIndexLocal(dm,dim,nCol,posCol,ic);
176:   MatSetValuesLocal(mat,nRow,ir,nCol,ic,val,insertMode);
177:   PetscFree2(ir,ic);
178:   return(0);
179: }

181: /*@C
182:   DMStagStencilToIndexLocal - Convert an array of DMStagStencil objects to an array of indices into a local vector.

184:   Not Collective

186:   Input Parameters:
187: + dm - the DMStag object
188: . dim - the dimension of the DMStag object
189: . n - the number of DMStagStencil objects
190: - pos - an array of n DMStagStencil objects

192:   Output Parameter:
193: . ix - output array of n indices

195:   Notes:
196:   The DMStagStencil objects in pos use global element indices.

198:   The .c fields in pos must always be set (even if to 0).

200:   Developer Notes:
201:   This is a "hot" function, and accepts the dimension redundantly to avoid having to perform any error checking inside the function.

203:   Level: developer

205: .seealso: DMSTAG, DMStagStencilLocation, DMStagStencil, DMGetLocalVector, DMCreateLocalVector
206: @*/
207: PetscErrorCode DMStagStencilToIndexLocal(DM dm,PetscInt dim,PetscInt n,const DMStagStencil *pos,PetscInt *ix)
208: {
209:   const DM_Stag * const stag = (DM_Stag*)dm->data;
210:   const PetscInt        epe = stag->entriesPerElement;

213:   if (dim == 1) {
214:     for (PetscInt idx=0; idx<n; ++idx) {
215:       const PetscInt eLocal = pos[idx].i - stag->startGhost[0];

217:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
218:     }
219:   } else if (dim == 2) {
220:     const PetscInt epr = stag->nGhost[0];

222:     for (PetscInt idx=0; idx<n; ++idx) {
223:       const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
224:       const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
225:       const PetscInt eLocal = eLocalx + epr*eLocaly;

227:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
228:     }
229:   } else if (dim == 3) {
230:     const PetscInt epr = stag->nGhost[0];
231:     const PetscInt epl = stag->nGhost[0]*stag->nGhost[1];

233:     for (PetscInt idx=0; idx<n; ++idx) {
234:       const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
235:       const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
236:       const PetscInt eLocalz = pos[idx].k - stag->startGhost[2];
237:       const PetscInt eLocal  = epl*eLocalz + epr*eLocaly + eLocalx;

239:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
240:     }
241:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
242:   return(0);
243: }

245: /*@C
246:   DMStagVecGetValuesStencil - get vector values using grid indexing

248:   Not Collective

250:   Input Parameters:
251: + dm - the DMStag object
252: . vec - the vector object
253: . n - the number of values to obtain
254: - pos - locations to obtain values from (as an array of DMStagStencil values)

256:   Output Parameter:
257: . val - value at the point

259:   Notes:
260:   Accepts stencils which refer to global element numbers, but
261:   only allows access to entries in the local representation (including ghosts).

263:   This approach is not as efficient as getting values directly with DMStagVecGetArray(), which is recommended for matrix free operators.

265:   Level: advanced

267: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), DMStagVecGetArray()
268: @*/
269: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec,PetscInt n,const DMStagStencil *pos,PetscScalar *val)
270: {
271:   PetscErrorCode    ierr;
272:   DM_Stag * const   stag = (DM_Stag*)dm->data;
273:   PetscInt          nLocal,dim,idx;
274:   PetscInt          *ix;
275:   PetscScalar const *arr;

280:   DMGetDimension(dm,&dim);
281:   VecGetLocalSize(vec,&nLocal);
282:   if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Vector should be a local vector. Local size %d does not match expected %d\n",nLocal,stag->entriesGhost);
283:   PetscMalloc1(n,&ix);
284:   DMStagStencilToIndexLocal(dm,dim,n,pos,ix);
285:   VecGetArrayRead(vec,&arr);
286:   for (idx=0; idx<n; ++idx) val[idx] = arr[ix[idx]];
287:   VecRestoreArrayRead(vec,&arr);
288:   PetscFree(ix);
289:   return(0);
290: }

292: /*@C
293:   DMStagVecSetValuesStencil - Set Vec values using global grid indexing

295:   Not Collective

297:   Input Parameters:
298: + dm - the DMStag object
299: . vec - the Vec
300: . n - the number of values to set
301: . pos - the locations to set values, as an array of DMStagStencil structs
302: . val - the values to set
303: - insertMode - INSERT_VALUES or ADD_VALUES

305:   Notes:
306:   The vector is expected to be a global vector compatible with the DM (usually obtained by DMGetGlobalVector() or DMCreateGlobalVector()).

308:   This approach is not as efficient as setting values directly with DMStagVecGetArray(), which is recommended for matrix-free operators.
309:   For assembling systems, where overhead may be less important than convenience, this routine could be helpful in assembling a righthand side and a matrix (using DMStagMatSetValuesStencil()).

311:   Level: advanced

313: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagMatSetValuesStencil(), DMCreateGlobalVector(), DMGetLocalVector(), DMStagVecGetArray()
314: @*/
315: PetscErrorCode DMStagVecSetValuesStencil(DM dm,Vec vec,PetscInt n,const DMStagStencil *pos,const PetscScalar *val,InsertMode insertMode)
316: {
317:   PetscErrorCode  ierr;
318:   DM_Stag * const stag = (DM_Stag*)dm->data;
319:   PetscInt        dim,nLocal;
320:   PetscInt        *ix;

325:   DMGetDimension(dm,&dim);
326:   VecGetLocalSize(vec,&nLocal);
327:   if (nLocal != stag->entries) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONG,"Provided vec has a different number of local entries (%D) than expected (%D). It should be a global vector",nLocal,stag->entries);
328:   PetscMalloc1(n,&ix);
329:   DMStagStencilToIndexLocal(dm,dim,n,pos,ix);
330:   VecSetValuesLocal(vec,n,ix,val,insertMode);
331:   PetscFree(ix);
332:   return(0);
333: }