Actual source code: stagstencil.c
petsc-3.13.6 2020-09-29
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: DMStagMatGetValuesStencil - retrieve local 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
176: Output Parameter:
177: . val - logically two-dimensional array of values
179: Level: advanced
181: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
182: @*/
183: PetscErrorCode DMStagMatGetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,PetscScalar *val)
184: {
186: PetscInt dim;
187: PetscInt *ir,*ic;
192: DMGetDimension(dm,&dim);
193: PetscMalloc2(nRow,&ir,nCol,&ic);
194: DMStagStencilToIndexLocal(dm,nRow,posRow,ir);
195: DMStagStencilToIndexLocal(dm,nCol,posCol,ic);
196: MatGetValuesLocal(mat,nRow,ir,nCol,ic,val);
197: PetscFree2(ir,ic);
198: return(0);
199: }
201: /*@C
202: DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing
204: Not Collective
206: Input Parameters:
207: + dm - the DMStag object
208: . mat - the matrix
209: . nRow - number of rows
210: . posRow - grid locations (including components) of rows
211: . nCol - number of columns
212: . posCol - grid locations (including components) of columns
213: . val - logically two-dimensional array of values
214: - insertmode - INSERT_VALUES or ADD_VALUES
216: Notes:
217: See notes for MatSetValuesStencil()
219: Level: intermediate
221: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatGetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
222: @*/
223: PetscErrorCode DMStagMatSetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,const PetscScalar *val,InsertMode insertMode)
224: {
226: PetscInt dim;
227: PetscInt *ir,*ic;
232: DMGetDimension(dm,&dim);
233: PetscMalloc2(nRow,&ir,nCol,&ic);
234: DMStagStencilToIndexLocal(dm,nRow,posRow,ir);
235: DMStagStencilToIndexLocal(dm,nCol,posCol,ic);
236: MatSetValuesLocal(mat,nRow,ir,nCol,ic,val,insertMode);
237: PetscFree2(ir,ic);
238: return(0);
239: }
241: /*@C
242: DMStagVecGetValuesStencil - get vector values using grid indexing
244: Not Collective
246: Input Parameters:
247: + dm - the DMStag object
248: . vec - the vector object
249: . n - the number of values to obtain
250: - pos - locations to obtain values from (as an array of DMStagStencil values)
252: Output Parameter:
253: . val - value at the point
255: Notes:
256: Accepts stencils which refer to global element numbers, but
257: only allows access to entries in the local representation (including ghosts).
259: This approach is not as efficient as setting values directly with DMStagVecGetArray(), which is recommended for matrix free operators.
261: Level: advanced
263: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), DMStagVecGetArray()
264: @*/
265: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec,PetscInt n,const DMStagStencil *pos,PetscScalar *val)
266: {
267: PetscErrorCode ierr;
268: DM_Stag * const stag = (DM_Stag*)dm->data;
269: PetscInt nLocal,dim,idx;
270: PetscInt *ix;
271: PetscScalar const *arr;
276: DMGetDimension(dm,&dim);
277: VecGetLocalSize(vec,&nLocal);
278: 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);
279: PetscMalloc1(n,&ix);
280: DMStagStencilToIndexLocal(dm,n,pos,ix);
281: VecGetArrayRead(vec,&arr);
282: for (idx=0; idx<n; ++idx) val[idx] = arr[ix[idx]];
283: VecRestoreArrayRead(vec,&arr);
284: PetscFree(ix);
285: return(0);
286: }
288: /*@C
289: DMStagVecSetValuesStencil - Set Vec values using global grid indexing
291: Not Collective
293: Input Parameters:
294: + dm - the DMStag object
295: . vec - the Vec
296: . n - the number of values to set
297: . pos - the locations to set values, as an array of DMStagStencil structs
298: . val - the values to set
299: - insertMode - INSERT_VALUES or ADD_VALUES
301: Notes:
302: The vector is expected to be a global vector compatible with the DM (usually obtained by DMGetGlobalVector() or DMCreateGlobalVector()).
304: This approach is not as efficient as setting values directly with DMStagVecGetArray(), which is recommended for matrix-free operators.
305: 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()).
307: Level: advanced
309: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagMatSetValuesStencil(), DMCreateGlobalVector(), DMGetLocalVector(), DMStagVecGetArray()
310: @*/
311: PetscErrorCode DMStagVecSetValuesStencil(DM dm,Vec vec,PetscInt n,const DMStagStencil *pos,const PetscScalar *val,InsertMode insertMode)
312: {
313: PetscErrorCode ierr;
314: DM_Stag * const stag = (DM_Stag*)dm->data;
315: PetscInt dim,nLocal;
316: PetscInt *ix;
321: DMGetDimension(dm,&dim);
322: VecGetLocalSize(vec,&nLocal);
323: 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);
324: PetscMalloc1(n,&ix);
325: DMStagStencilToIndexLocal(dm,n,pos,ix);
326: VecSetValuesLocal(vec,n,ix,val,insertMode);
327: PetscFree(ix);
328: return(0);
329: }