Actual source code: stagstencil.c
petsc-3.14.6 2021-03-30
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: }