Actual source code: stagstencil.c

petsc-master 2019-11-11
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 defined(PETSC_USE_DEBUG)
116:   {
117:     PetscInt i,nGhost[DMSTAG_MAX_DIM],endGhost[DMSTAG_MAX_DIM];
118:     DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],&nGhost[0],&nGhost[1],&nGhost[2]);
119:     for (i=0; i<DMSTAG_MAX_DIM; ++i) endGhost[i] = startGhost[i] + nGhost[i];
120:     for (i=0; i<n; ++i) {
121:       PetscInt dof;
122:       DMStagGetLocationDOF(dm,pos[i].loc,&dof);
123:       if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[pos[i].loc]);
124:       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);
125:       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);
126:       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);
127:       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);
128:       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);
129:     }
130:   }
131: #else
132:   DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
133: #endif
134:   if (dim == 1) {
135:     for (idx=0; idx<n; ++idx) {
136:       const PetscInt eLocal = pos[idx].i - startGhost[0]; /* Local element number */
137:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
138:     }
139:   } else if (dim == 2) {
140:     const PetscInt epr = stag->nGhost[0];
141:     DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],NULL,NULL,NULL,NULL);
142:     for (idx=0; idx<n; ++idx) {
143:       const PetscInt eLocalx = pos[idx].i - startGhost[0];
144:       const PetscInt eLocaly = pos[idx].j - startGhost[1];
145:       const PetscInt eLocal = eLocalx + epr*eLocaly;
146:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
147:     }
148:   } else if (dim == 3) {
149:     const PetscInt epr = stag->nGhost[0];
150:     const PetscInt epl = stag->nGhost[0]*stag->nGhost[1];
151:     DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);
152:     for (idx=0; idx<n; ++idx) {
153:       const PetscInt eLocalx = pos[idx].i - startGhost[0];
154:       const PetscInt eLocaly = pos[idx].j - startGhost[1];
155:       const PetscInt eLocalz = pos[idx].k - startGhost[2];
156:       const PetscInt eLocal  = epl*eLocalz + epr*eLocaly + eLocalx;
157:       ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
158:     }
159:   } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
160:   return(0);
161: }

163: /*@C
164:   DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing

166:   Not Collective

168:   Input Parameters:
169: + dm - the DMStag object
170: . mat - the matrix
171: . nRow - number of rows
172: . posRow - grid locations (including components) of rows
173: . nCol - number of columns
174: . posCol - grid locations (including components) of columns
175: . val - logically two-dimensional array of values
176: - insertmode - INSERT_VALUES or ADD_VALUES

178:   Notes:
179:   See notes for MatSetValuesStencil()

181:   Level: intermediate

183: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
184: @*/
185: PetscErrorCode DMStagMatSetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,const PetscScalar *val,InsertMode insertMode)
186: {
188:   PetscInt       dim;
189:   PetscInt       *ir,*ic;

194:   DMGetDimension(dm,&dim);
195:   PetscMalloc2(nRow,&ir,nCol,&ic);
196:   DMStagStencilToIndexLocal(dm,nRow,posRow,ir);
197:   DMStagStencilToIndexLocal(dm,nCol,posCol,ic);
198:   MatSetValuesLocal(mat,nRow,ir,nCol,ic,val,insertMode);
199:   PetscFree2(ir,ic);
200:   return(0);
201: }

203: /*@C
204:   DMStagVecGetValuesStencil - get vector values using grid indexing

206:   Not Collective

208:   Input Parameters:
209: + dm - the DMStag object
210: . vec - the vector object
211: . n - the number of values to obtain
212: - pos - locations to obtain values from (as an array of DMStagStencil values)

214:   Output Parameter:
215: . val - value at the point

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

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

223:   Level: advanced

225: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), DMStagVecGetArrayDOF()
226: @*/
227: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec,PetscInt n,const DMStagStencil *pos,PetscScalar *val)
228: {
229:   PetscErrorCode    ierr;
230:   DM_Stag * const   stag = (DM_Stag*)dm->data;
231:   PetscInt          nLocal,dim,idx;
232:   PetscInt          *ix;
233:   PetscScalar const *arr;

238:   DMGetDimension(dm,&dim);
239:   VecGetLocalSize(vec,&nLocal);
240:   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);
241:   PetscMalloc1(n,&ix);
242:   DMStagStencilToIndexLocal(dm,n,pos,ix);
243:   VecGetArrayRead(vec,&arr);
244:   for (idx=0; idx<n; ++idx) val[idx] = arr[ix[idx]];
245:   VecRestoreArrayRead(vec,&arr);
246:   PetscFree(ix);
247:   return(0);
248: }

250: /*@C
251:   DMStagVecSetValuesStencil - Set Vec values using global grid indexing

253:   Not Collective

255:   Input Parameters:
256: + dm - the DMStag object
257: . vec - the Vec
258: . n - the number of values to set
259: . pos - the locations to set values, as an array of DMStagStencil structs
260: . val - the values to set
261: - insertMode - INSERT_VALUES or ADD_VALUES

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

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

269:   Level: advanced

271: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagMatSetValuesStencil(), DMCreateGlobalVector(), DMGetLocalVector(), DMStagVecGetArrayDOF()
272: @*/
273: PetscErrorCode DMStagVecSetValuesStencil(DM dm,Vec vec,PetscInt n,const DMStagStencil *pos,const PetscScalar *val,InsertMode insertMode)
274: {
275:   PetscErrorCode  ierr;
276:   DM_Stag * const stag = (DM_Stag*)dm->data;
277:   PetscInt        dim,nLocal;
278:   PetscInt        *ix;

283:   DMGetDimension(dm,&dim);
284:   VecGetLocalSize(vec,&nLocal);
285:   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);
286:   PetscMalloc1(n,&ix);
287:   DMStagStencilToIndexLocal(dm,n,pos,ix);
288:   VecSetValuesLocal(vec,n,ix,val,insertMode);
289:   PetscFree(ix);
290:   return(0);
291: }