Actual source code: stagstencil.c
petsc-3.11.4 2019-09-28
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
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: }