Actual source code: stagstencil.c

petsc-master 2020-08-10
Report Typos and Errors
  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"};
  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: /* Convert an array of DMStagStencil objects to an array of indices into a local vector.
104:   The .c fields in pos must always be set (even if to 0).  */
105: static PetscErrorCode DMStagStencilToIndexLocal(DM dm,PetscInt n,const DMStagStencil *pos,PetscInt *ix)
106: {
107:   PetscErrorCode        ierr;
108:   const DM_Stag * const stag = (DM_Stag*)dm->data;
109:   PetscInt              idx,dim,startGhost[DMSTAG_MAX_DIM];
110:   const PetscInt        epe = stag->entriesPerElement;

114:   DMGetDimension(dm,&dim);
115:   if (PetscDefined(USE_DEBUG)) {
116:     PetscInt i,nGhost[DMSTAG_MAX_DIM],endGhost[DMSTAG_MAX_DIM];
117:     DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],&nGhost[0],&nGhost[1],&nGhost[2]);
118:     for (i=0; i<DMSTAG_MAX_DIM; ++i) endGhost[i] = startGhost[i] + nGhost[i];
119:     for (i=0; i<n; ++i) {
120:       PetscInt dof;
121:       DMStagGetLocationDOF(dm,pos[i].loc,&dof);
122:       if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[pos[i].loc]);
123:       if (pos[i].c < 0) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Negative component number (%d) supplied in loc[%D]",pos[i].c,i);
124:       if (pos[i].c > dof-1) SETERRQ3(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Supplied component number (%D) for location %s is too big (maximum %D)",pos[i].c,DMStagStencilLocations[pos[i].loc],dof-1);
125:       if (            pos[i].i >= endGhost[0] || pos[i].i < startGhost[0] ) SETERRQ3(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Supplied x element index %D out of range. Should be in [%D,%D]",pos[i].i,startGhost[0],endGhost[0]-1);
126:       if (dim > 1 && (pos[i].j >= endGhost[1] || pos[i].j < startGhost[1])) SETERRQ3(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Supplied y element index %D out of range. Should be in [%D,%D]",pos[i].j,startGhost[1],endGhost[1]-1);
127:       if (dim > 2 && (pos[i].k >= endGhost[2] || pos[i].k < startGhost[2])) SETERRQ3(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Supplied z element index %D out of range. Should be in [%D,%D]",pos[i].k,startGhost[2],endGhost[2]-1);
128:     }
129:   } else {
130:     DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
131:   }
132:   if (dim == 1) {
133:     for (idx=0; idx<n; ++idx) {
134:       const PetscInt eLocal = pos[idx].i - startGhost[0]; /* Local element number */
135:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
136:     }
137:   } else if (dim == 2) {
138:     const PetscInt epr = stag->nGhost[0];
139:     DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],NULL,NULL,NULL,NULL);
140:     for (idx=0; idx<n; ++idx) {
141:       const PetscInt eLocalx = pos[idx].i - startGhost[0];
142:       const PetscInt eLocaly = pos[idx].j - startGhost[1];
143:       const PetscInt eLocal = eLocalx + epr*eLocaly;
144:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
145:     }
146:   } else if (dim == 3) {
147:     const PetscInt epr = stag->nGhost[0];
148:     const PetscInt epl = stag->nGhost[0]*stag->nGhost[1];
149:     DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
150:     for (idx=0; idx<n; ++idx) {
151:       const PetscInt eLocalx = pos[idx].i - startGhost[0];
152:       const PetscInt eLocaly = pos[idx].j - startGhost[1];
153:       const PetscInt eLocalz = pos[idx].k - startGhost[2];
154:       const PetscInt eLocal  = epl*eLocalz + epr*eLocaly + eLocalx;
155:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
156:     }
157:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
158:   return(0);
159: }

161: /*@C
162:   DMStagMatGetValuesStencil - retrieve local matrix entries using grid indexing

164:   Not Collective

166:   Input Parameters:
167: + dm - the DMStag object
168: . mat - the matrix
169: . nRow - number of rows
170: . posRow - grid locations (including components) of rows
171: . nCol - number of columns
172: - posCol - grid locations (including components) of columns

174:   Output Parameter:
175: . val - logically two-dimensional array of values

177:   Level: advanced

179: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
180: @*/
181: PetscErrorCode DMStagMatGetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,PetscScalar *val)
182: {
184:   PetscInt       dim;
185:   PetscInt       *ir,*ic;

190:   DMGetDimension(dm,&dim);
191:   PetscMalloc2(nRow,&ir,nCol,&ic);
192:   DMStagStencilToIndexLocal(dm,nRow,posRow,ir);
193:   DMStagStencilToIndexLocal(dm,nCol,posCol,ic);
194:   MatGetValuesLocal(mat,nRow,ir,nCol,ic,val);
195:   PetscFree2(ir,ic);
196:   return(0);
197: }

199: /*@C
200:   DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing

202:   Not Collective

204:   Input Parameters:
205: + dm - the DMStag object
206: . mat - the matrix
207: . nRow - number of rows
208: . posRow - grid locations (including components) of rows
209: . nCol - number of columns
210: . posCol - grid locations (including components) of columns
211: . val - logically two-dimensional array of values
212: - insertmode - INSERT_VALUES or ADD_VALUES

214:   Notes:
215:   See notes for MatSetValuesStencil()

217:   Level: intermediate

219: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatGetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
220: @*/
221: PetscErrorCode DMStagMatSetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,const PetscScalar *val,InsertMode insertMode)
222: {
224:   PetscInt       dim;
225:   PetscInt       *ir,*ic;

230:   DMGetDimension(dm,&dim);
231:   PetscMalloc2(nRow,&ir,nCol,&ic);
232:   DMStagStencilToIndexLocal(dm,nRow,posRow,ir);
233:   DMStagStencilToIndexLocal(dm,nCol,posCol,ic);
234:   MatSetValuesLocal(mat,nRow,ir,nCol,ic,val,insertMode);
235:   PetscFree2(ir,ic);
236:   return(0);
237: }

239: /*@C
240:   DMStagVecGetValuesStencil - get vector values using grid indexing

242:   Not Collective

244:   Input Parameters:
245: + dm - the DMStag object
246: . vec - the vector object
247: . n - the number of values to obtain
248: - pos - locations to obtain values from (as an array of DMStagStencil values)

250:   Output Parameter:
251: . val - value at the point

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

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

259:   Level: advanced

261: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), DMStagVecGetArray()
262: @*/
263: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec,PetscInt n,const DMStagStencil *pos,PetscScalar *val)
264: {
265:   PetscErrorCode    ierr;
266:   DM_Stag * const   stag = (DM_Stag*)dm->data;
267:   PetscInt          nLocal,dim,idx;
268:   PetscInt          *ix;
269:   PetscScalar const *arr;

274:   DMGetDimension(dm,&dim);
275:   VecGetLocalSize(vec,&nLocal);
276:   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);
277:   PetscMalloc1(n,&ix);
278:   DMStagStencilToIndexLocal(dm,n,pos,ix);
279:   VecGetArrayRead(vec,&arr);
280:   for (idx=0; idx<n; ++idx) val[idx] = arr[ix[idx]];
281:   VecRestoreArrayRead(vec,&arr);
282:   PetscFree(ix);
283:   return(0);
284: }

286: /*@C
287:   DMStagVecSetValuesStencil - Set Vec values using global grid indexing

289:   Not Collective

291:   Input Parameters:
292: + dm - the DMStag object
293: . vec - the Vec
294: . n - the number of values to set
295: . pos - the locations to set values, as an array of DMStagStencil structs
296: . val - the values to set
297: - insertMode - INSERT_VALUES or ADD_VALUES

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

302:   This approach is not as efficient as setting values directly with DMStagVecGetArray(), which is recommended for matrix-free operators. 
303:   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()).

305:   Level: advanced

307: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagMatSetValuesStencil(), DMCreateGlobalVector(), DMGetLocalVector(), DMStagVecGetArray()
308: @*/
309: PetscErrorCode DMStagVecSetValuesStencil(DM dm,Vec vec,PetscInt n,const DMStagStencil *pos,const PetscScalar *val,InsertMode insertMode)
310: {
311:   PetscErrorCode  ierr;
312:   DM_Stag * const stag = (DM_Stag*)dm->data;
313:   PetscInt        dim,nLocal;
314:   PetscInt        *ix;

319:   DMGetDimension(dm,&dim);
320:   VecGetLocalSize(vec,&nLocal);
321:   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);
322:   PetscMalloc1(n,&ix);
323:   DMStagStencilToIndexLocal(dm,n,pos,ix);
324:   VecSetValuesLocal(vec,n,ix,val,insertMode);
325:   PetscFree(ix);
326:   return(0);
327: }