Actual source code: stagstencil.c
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: /*@C
104: DMStagMatGetValuesStencil - retrieve local matrix entries using grid indexing
106: Not Collective
108: Input Parameters:
109: + dm - the DMStag object
110: . mat - the matrix
111: . nRow - number of rows
112: . posRow - grid locations (including components) of rows
113: . nCol - number of columns
114: - posCol - grid locations (including components) of columns
116: Output Parameter:
117: . val - logically two-dimensional array of values
119: Level: advanced
121: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
122: @*/
123: PetscErrorCode DMStagMatGetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,PetscScalar *val)
124: {
126: PetscInt dim;
127: PetscInt *ir,*ic;
132: DMGetDimension(dm,&dim);
133: PetscMalloc2(nRow,&ir,nCol,&ic);
134: DMStagStencilToIndexLocal(dm,dim,nRow,posRow,ir);
135: DMStagStencilToIndexLocal(dm,dim,nCol,posCol,ic);
136: MatGetValuesLocal(mat,nRow,ir,nCol,ic,val);
137: PetscFree2(ir,ic);
138: return(0);
139: }
141: /*@C
142: DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing
144: Not Collective
146: Input Parameters:
147: + dm - the DMStag object
148: . mat - the matrix
149: . nRow - number of rows
150: . posRow - grid locations (including components) of rows
151: . nCol - number of columns
152: . posCol - grid locations (including components) of columns
153: . val - logically two-dimensional array of values
154: - insertmode - INSERT_VALUES or ADD_VALUES
156: Notes:
157: See notes for MatSetValuesStencil()
159: Level: intermediate
161: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatGetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
162: @*/
163: PetscErrorCode DMStagMatSetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,const PetscScalar *val,InsertMode insertMode)
164: {
166: PetscInt dim;
167: PetscInt *ir,*ic;
172: DMGetDimension(dm,&dim);
173: PetscMalloc2(nRow,&ir,nCol,&ic);
174: DMStagStencilToIndexLocal(dm,dim,nRow,posRow,ir);
175: DMStagStencilToIndexLocal(dm,dim,nCol,posCol,ic);
176: MatSetValuesLocal(mat,nRow,ir,nCol,ic,val,insertMode);
177: PetscFree2(ir,ic);
178: return(0);
179: }
181: /*@C
182: DMStagStencilToIndexLocal - Convert an array of DMStagStencil objects to an array of indices into a local vector.
184: Not Collective
186: Input Parameters:
187: + dm - the DMStag object
188: . dim - the dimension of the DMStag object
189: . n - the number of DMStagStencil objects
190: - pos - an array of n DMStagStencil objects
192: Output Parameter:
193: . ix - output array of n indices
195: Notes:
196: The DMStagStencil objects in pos use global element indices.
198: The .c fields in pos must always be set (even if to 0).
200: Developer Notes:
201: This is a "hot" function, and accepts the dimension redundantly to avoid having to perform any error checking inside the function.
203: Level: developer
205: .seealso: DMSTAG, DMStagStencilLocation, DMStagStencil, DMGetLocalVector, DMCreateLocalVector
206: @*/
207: PetscErrorCode DMStagStencilToIndexLocal(DM dm,PetscInt dim,PetscInt n,const DMStagStencil *pos,PetscInt *ix)
208: {
209: const DM_Stag * const stag = (DM_Stag*)dm->data;
210: const PetscInt epe = stag->entriesPerElement;
213: if (dim == 1) {
214: for (PetscInt idx=0; idx<n; ++idx) {
215: const PetscInt eLocal = pos[idx].i - stag->startGhost[0];
217: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
218: }
219: } else if (dim == 2) {
220: const PetscInt epr = stag->nGhost[0];
222: for (PetscInt idx=0; idx<n; ++idx) {
223: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
224: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
225: const PetscInt eLocal = eLocalx + epr*eLocaly;
227: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
228: }
229: } else if (dim == 3) {
230: const PetscInt epr = stag->nGhost[0];
231: const PetscInt epl = stag->nGhost[0]*stag->nGhost[1];
233: for (PetscInt idx=0; idx<n; ++idx) {
234: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
235: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
236: const PetscInt eLocalz = pos[idx].k - stag->startGhost[2];
237: const PetscInt eLocal = epl*eLocalz + epr*eLocaly + eLocalx;
239: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
240: }
241: } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
242: return(0);
243: }
245: /*@C
246: DMStagVecGetValuesStencil - get vector values using grid indexing
248: Not Collective
250: Input Parameters:
251: + dm - the DMStag object
252: . vec - the vector object
253: . n - the number of values to obtain
254: - pos - locations to obtain values from (as an array of DMStagStencil values)
256: Output Parameter:
257: . val - value at the point
259: Notes:
260: Accepts stencils which refer to global element numbers, but
261: only allows access to entries in the local representation (including ghosts).
263: This approach is not as efficient as setting values directly with DMStagVecGetArray(), which is recommended for matrix free operators.
265: Level: advanced
267: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), DMStagVecGetArray()
268: @*/
269: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec,PetscInt n,const DMStagStencil *pos,PetscScalar *val)
270: {
271: PetscErrorCode ierr;
272: DM_Stag * const stag = (DM_Stag*)dm->data;
273: PetscInt nLocal,dim,idx;
274: PetscInt *ix;
275: PetscScalar const *arr;
280: DMGetDimension(dm,&dim);
281: VecGetLocalSize(vec,&nLocal);
282: 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);
283: PetscMalloc1(n,&ix);
284: DMStagStencilToIndexLocal(dm,dim,n,pos,ix);
285: VecGetArrayRead(vec,&arr);
286: for (idx=0; idx<n; ++idx) val[idx] = arr[ix[idx]];
287: VecRestoreArrayRead(vec,&arr);
288: PetscFree(ix);
289: return(0);
290: }
292: /*@C
293: DMStagVecSetValuesStencil - Set Vec values using global grid indexing
295: Not Collective
297: Input Parameters:
298: + dm - the DMStag object
299: . vec - the Vec
300: . n - the number of values to set
301: . pos - the locations to set values, as an array of DMStagStencil structs
302: . val - the values to set
303: - insertMode - INSERT_VALUES or ADD_VALUES
305: Notes:
306: The vector is expected to be a global vector compatible with the DM (usually obtained by DMGetGlobalVector() or DMCreateGlobalVector()).
308: This approach is not as efficient as setting values directly with DMStagVecGetArray(), which is recommended for matrix-free operators.
309: 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()).
311: Level: advanced
313: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagMatSetValuesStencil(), DMCreateGlobalVector(), DMGetLocalVector(), DMStagVecGetArray()
314: @*/
315: PetscErrorCode DMStagVecSetValuesStencil(DM dm,Vec vec,PetscInt n,const DMStagStencil *pos,const PetscScalar *val,InsertMode insertMode)
316: {
317: PetscErrorCode ierr;
318: DM_Stag * const stag = (DM_Stag*)dm->data;
319: PetscInt dim,nLocal;
320: PetscInt *ix;
325: DMGetDimension(dm,&dim);
326: VecGetLocalSize(vec,&nLocal);
327: 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);
328: PetscMalloc1(n,&ix);
329: DMStagStencilToIndexLocal(dm,dim,n,pos,ix);
330: VecSetValuesLocal(vec,n,ix,val,insertMode);
331: PetscFree(ix);
332: return(0);
333: }