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", "DMStagStencilLocation", "", NULL};
10: /*@C
11: DMStagCreateISFromStencils - Create an `IS`, using global numberings, for a subset of DOF in a `DMSTAG` object
13: Collective
15: Input Parameters:
16: + dm - the `DMStag` object
17: . n_stencil - the number of stencils provided
18: - stencils - an array of `DMStagStencil` objects (`i`, `j`, and `k` are ignored)
20: Output Parameter:
21: . is - the global IS
23: Note:
24: Redundant entries in the stencils argument are ignored
26: Level: advanced
28: .seealso: [](chapter_stag), `DMSTAG`, `IS`, `DMStagStencil`, `DMCreateGlobalVector`
29: @*/
30: PetscErrorCode DMStagCreateISFromStencils(DM dm, PetscInt n_stencil, DMStagStencil *stencils, IS *is)
31: {
32: PetscInt *stencil_active;
33: DMStagStencil *stencils_ordered_unique;
34: PetscInt *idx, *idxLocal;
35: const PetscInt *ltogidx;
36: PetscInt n_stencil_unique, dim, count, nidx, nc_max;
37: ISLocalToGlobalMapping ltog;
38: PetscInt start[DMSTAG_MAX_DIM], n[DMSTAG_MAX_DIM], extraPoint[DMSTAG_MAX_DIM];
40: DMGetDimension(dm, &dim);
43: /* To assert that the resulting IS has unique, sorted, entries, we perform
44: a bucket sort, taking advantage of the fact that DMStagStencilLocation
45: enum values are integers starting with 1, in canonical order */
46: nc_max = 1; // maximum number of components to represent these stencils
47: n_stencil_unique = 0;
48: for (PetscInt p = 0; p < n_stencil; ++p) nc_max = PetscMax(nc_max, (stencils[p].c + 1));
49: PetscCalloc1(DMSTAG_NUMBER_LOCATIONS * nc_max, &stencil_active);
50: for (PetscInt p = 0; p < n_stencil; ++p) {
51: DMStagStencilLocation loc_canonical;
52: PetscInt slot;
54: DMStagStencilLocationCanonicalize(stencils[p].loc, &loc_canonical);
55: slot = nc_max * ((PetscInt)loc_canonical) + stencils[p].c;
56: if (stencil_active[slot] == 0) {
57: stencil_active[slot] = 1;
58: ++n_stencil_unique;
59: }
60: }
61: PetscMalloc1(n_stencil_unique, &stencils_ordered_unique);
62: {
63: PetscInt p = 0;
65: for (PetscInt i = 1; i < DMSTAG_NUMBER_LOCATIONS; ++i) {
66: for (PetscInt c = 0; c < nc_max; ++c) {
67: if (stencil_active[nc_max * i + c] != 0) {
68: stencils_ordered_unique[p].loc = (DMStagStencilLocation)i;
69: stencils_ordered_unique[p].c = c;
70: ++p;
71: }
72: }
73: }
74: }
75: PetscFree(stencil_active);
77: PetscMalloc1(n_stencil_unique, &idxLocal);
78: DMGetLocalToGlobalMapping(dm, <og);
79: ISLocalToGlobalMappingGetIndices(ltog, <ogidx);
80: DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &extraPoint[0], &extraPoint[1], &extraPoint[2]);
81: for (PetscInt d = dim; d < DMSTAG_MAX_DIM; ++d) {
82: start[d] = 0;
83: n[d] = 1; /* To allow for a single loop nest below */
84: extraPoint[d] = 0;
85: }
86: nidx = n_stencil_unique;
87: for (PetscInt d = 0; d < dim; ++d) nidx *= (n[d] + 1); /* Overestimate (always assumes extraPoint) */
88: PetscMalloc1(nidx, &idx);
89: count = 0;
90: /* Note that unused loop variables are not accessed, for lower dimensions */
91: for (PetscInt k = start[2]; k < start[2] + n[2] + extraPoint[2]; ++k) {
92: for (PetscInt j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
93: for (PetscInt i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
94: for (PetscInt p = 0; p < n_stencil_unique; ++p) {
95: stencils_ordered_unique[p].i = i;
96: stencils_ordered_unique[p].j = j;
97: stencils_ordered_unique[p].k = k;
98: }
99: DMStagStencilToIndexLocal(dm, dim, n_stencil_unique, stencils_ordered_unique, idxLocal);
100: for (PetscInt p = 0; p < n_stencil_unique; ++p) {
101: const PetscInt gidx = ltogidx[idxLocal[p]];
102: if (gidx >= 0) {
103: idx[count] = gidx;
104: ++count;
105: }
106: }
107: }
108: }
109: }
111: ISLocalToGlobalMappingRestoreIndices(ltog, <ogidx);
112: PetscFree(stencils_ordered_unique);
113: PetscFree(idxLocal);
115: ISCreateGeneral(PetscObjectComm((PetscObject)dm), count, idx, PETSC_OWN_POINTER, is);
116: ISSetInfo(*is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE);
117: ISSetInfo(*is, IS_UNIQUE, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE);
119: return 0;
120: }
122: /*@C
123: DMStagGetLocationDOF - Get number of DOF associated with a given point in a `DMSTAG` grid
125: Not Collective
127: Input Parameters:
128: + dm - the `DMSTAG` object
129: - loc - grid point (see `DMStagStencilLocation`)
131: Output Parameter:
132: . dof - the number of DOF (components) living at `loc` in `dm`
134: Level: intermediate
136: .seealso: [](chapter_stag), `DMSTAG`, `DMStagStencilLocation`, `DMStagStencil`, `DMDAGetDof()`
137: @*/
138: PetscErrorCode DMStagGetLocationDOF(DM dm, DMStagStencilLocation loc, PetscInt *dof)
139: {
140: const DM_Stag *const stag = (DM_Stag *)dm->data;
141: PetscInt dim;
144: DMGetDimension(dm, &dim);
145: switch (dim) {
146: case 1:
147: switch (loc) {
148: case DMSTAG_LEFT:
149: case DMSTAG_RIGHT:
150: *dof = stag->dof[0];
151: break;
152: case DMSTAG_ELEMENT:
153: *dof = stag->dof[1];
154: break;
155: default:
156: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
157: }
158: break;
159: case 2:
160: switch (loc) {
161: case DMSTAG_DOWN_LEFT:
162: case DMSTAG_DOWN_RIGHT:
163: case DMSTAG_UP_LEFT:
164: case DMSTAG_UP_RIGHT:
165: *dof = stag->dof[0];
166: break;
167: case DMSTAG_LEFT:
168: case DMSTAG_RIGHT:
169: case DMSTAG_UP:
170: case DMSTAG_DOWN:
171: *dof = stag->dof[1];
172: break;
173: case DMSTAG_ELEMENT:
174: *dof = stag->dof[2];
175: break;
176: default:
177: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
178: }
179: break;
180: case 3:
181: switch (loc) {
182: case DMSTAG_BACK_DOWN_LEFT:
183: case DMSTAG_BACK_DOWN_RIGHT:
184: case DMSTAG_BACK_UP_LEFT:
185: case DMSTAG_BACK_UP_RIGHT:
186: case DMSTAG_FRONT_DOWN_LEFT:
187: case DMSTAG_FRONT_DOWN_RIGHT:
188: case DMSTAG_FRONT_UP_LEFT:
189: case DMSTAG_FRONT_UP_RIGHT:
190: *dof = stag->dof[0];
191: break;
192: case DMSTAG_BACK_DOWN:
193: case DMSTAG_BACK_LEFT:
194: case DMSTAG_BACK_RIGHT:
195: case DMSTAG_BACK_UP:
196: case DMSTAG_DOWN_LEFT:
197: case DMSTAG_DOWN_RIGHT:
198: case DMSTAG_UP_LEFT:
199: case DMSTAG_UP_RIGHT:
200: case DMSTAG_FRONT_DOWN:
201: case DMSTAG_FRONT_LEFT:
202: case DMSTAG_FRONT_RIGHT:
203: case DMSTAG_FRONT_UP:
204: *dof = stag->dof[1];
205: break;
206: case DMSTAG_LEFT:
207: case DMSTAG_RIGHT:
208: case DMSTAG_DOWN:
209: case DMSTAG_UP:
210: case DMSTAG_BACK:
211: case DMSTAG_FRONT:
212: *dof = stag->dof[2];
213: break;
214: case DMSTAG_ELEMENT:
215: *dof = stag->dof[3];
216: break;
217: default:
218: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
219: }
220: break;
221: default:
222: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
223: }
224: return 0;
225: }
227: /*
228: Convert to a location value with only BACK, DOWN, LEFT, and ELEMENT involved
229: */
230: PETSC_INTERN PetscErrorCode DMStagStencilLocationCanonicalize(DMStagStencilLocation loc, DMStagStencilLocation *locCanonical)
231: {
232: switch (loc) {
233: case DMSTAG_ELEMENT:
234: *locCanonical = DMSTAG_ELEMENT;
235: break;
236: case DMSTAG_LEFT:
237: case DMSTAG_RIGHT:
238: *locCanonical = DMSTAG_LEFT;
239: break;
240: case DMSTAG_DOWN:
241: case DMSTAG_UP:
242: *locCanonical = DMSTAG_DOWN;
243: break;
244: case DMSTAG_BACK:
245: case DMSTAG_FRONT:
246: *locCanonical = DMSTAG_BACK;
247: break;
248: case DMSTAG_DOWN_LEFT:
249: case DMSTAG_DOWN_RIGHT:
250: case DMSTAG_UP_LEFT:
251: case DMSTAG_UP_RIGHT:
252: *locCanonical = DMSTAG_DOWN_LEFT;
253: break;
254: case DMSTAG_BACK_LEFT:
255: case DMSTAG_BACK_RIGHT:
256: case DMSTAG_FRONT_LEFT:
257: case DMSTAG_FRONT_RIGHT:
258: *locCanonical = DMSTAG_BACK_LEFT;
259: break;
260: case DMSTAG_BACK_DOWN:
261: case DMSTAG_BACK_UP:
262: case DMSTAG_FRONT_DOWN:
263: case DMSTAG_FRONT_UP:
264: *locCanonical = DMSTAG_BACK_DOWN;
265: break;
266: case DMSTAG_BACK_DOWN_LEFT:
267: case DMSTAG_BACK_DOWN_RIGHT:
268: case DMSTAG_BACK_UP_LEFT:
269: case DMSTAG_BACK_UP_RIGHT:
270: case DMSTAG_FRONT_DOWN_LEFT:
271: case DMSTAG_FRONT_DOWN_RIGHT:
272: case DMSTAG_FRONT_UP_LEFT:
273: case DMSTAG_FRONT_UP_RIGHT:
274: *locCanonical = DMSTAG_BACK_DOWN_LEFT;
275: break;
276: default:
277: *locCanonical = DMSTAG_NULL_LOCATION;
278: break;
279: }
280: return 0;
281: }
283: /*@C
284: DMStagMatGetValuesStencil - retrieve local matrix entries using grid indexing
286: Not Collective
288: Input Parameters:
289: + dm - the `DMSTAG` object
290: . mat - the matrix
291: . nRow - number of rows
292: . posRow - grid locations (including components) of rows
293: . nCol - number of columns
294: - posCol - grid locations (including components) of columns
296: Output Parameter:
297: . val - logically two-dimensional array of values
299: Level: advanced
301: .seealso: [](chapter_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagVecSetValuesStencil()`, `DMStagMatSetValuesStencil()`, `MatSetValuesStencil()`, `MatAssemblyBegin()`, `MatAssemblyEnd()`, `DMCreateMatrix()`
302: @*/
303: PetscErrorCode DMStagMatGetValuesStencil(DM dm, Mat mat, PetscInt nRow, const DMStagStencil *posRow, PetscInt nCol, const DMStagStencil *posCol, PetscScalar *val)
304: {
305: PetscInt dim;
306: PetscInt *ir, *ic;
310: DMGetDimension(dm, &dim);
311: PetscMalloc2(nRow, &ir, nCol, &ic);
312: DMStagStencilToIndexLocal(dm, dim, nRow, posRow, ir);
313: DMStagStencilToIndexLocal(dm, dim, nCol, posCol, ic);
314: MatGetValuesLocal(mat, nRow, ir, nCol, ic, val);
315: PetscFree2(ir, ic);
316: return 0;
317: }
319: /*@C
320: DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing
322: Not Collective
324: Input Parameters:
325: + dm - the `DMSTAG` object
326: . mat - the matrix
327: . nRow - number of rows
328: . posRow - grid locations (including components) of rows
329: . nCol - number of columns
330: . posCol - grid locations (including components) of columns
331: . val - logically two-dimensional array of values
332: - insertmode - `INSERT_VALUES` or `ADD_VALUES`
334: Notes:
335: See notes for `MatSetValuesStencil()`
337: Level: intermediate
339: .seealso: [](chapter_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagVecSetValuesStencil()`, `DMStagMatGetValuesStencil()`, `MatSetValuesStencil()`, `MatAssemblyBegin()`, `MatAssemblyEnd()`, `DMCreateMatrix()`
340: @*/
341: PetscErrorCode DMStagMatSetValuesStencil(DM dm, Mat mat, PetscInt nRow, const DMStagStencil *posRow, PetscInt nCol, const DMStagStencil *posCol, const PetscScalar *val, InsertMode insertMode)
342: {
343: PetscInt *ir, *ic;
347: PetscMalloc2(nRow, &ir, nCol, &ic);
348: DMStagStencilToIndexLocal(dm, dm->dim, nRow, posRow, ir);
349: DMStagStencilToIndexLocal(dm, dm->dim, nCol, posCol, ic);
350: MatSetValuesLocal(mat, nRow, ir, nCol, ic, val, insertMode);
351: PetscFree2(ir, ic);
352: return 0;
353: }
355: /*@C
356: DMStagStencilToIndexLocal - Convert an array of `DMStagStenci`l objects to an array of indices into a local vector.
358: Not Collective
360: Input Parameters:
361: + dm - the `DMSTAG` object
362: . dim - the dimension of the `DMSTAG` object
363: . n - the number of `DMStagStencil` objects
364: - pos - an array of `n` `DMStagStencil` objects
366: Output Parameter:
367: . ix - output array of `n` indices
369: Notes:
370: The `DMStagStencil` objects in `pos` use global element indices.
372: The `.c` fields in `pos` must always be set (even if to `0`).
374: Developer Notes:
375: This is a "hot" function, and accepts the dimension redundantly to avoid having to perform any error checking inside the function.
377: Level: developer
379: .seealso: [](chapter_stag), `DMSTAG`, `DMStagStencilLocation`, `DMStagStencil`, `DMGetLocalVector`, `DMCreateLocalVector`
380: @*/
381: PetscErrorCode DMStagStencilToIndexLocal(DM dm, PetscInt dim, PetscInt n, const DMStagStencil *pos, PetscInt *ix)
382: {
383: const DM_Stag *const stag = (DM_Stag *)dm->data;
384: const PetscInt epe = stag->entriesPerElement;
387: if (dim == 1) {
388: for (PetscInt idx = 0; idx < n; ++idx) {
389: const PetscInt eLocal = pos[idx].i - stag->startGhost[0];
391: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
392: }
393: } else if (dim == 2) {
394: const PetscInt epr = stag->nGhost[0];
396: for (PetscInt idx = 0; idx < n; ++idx) {
397: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
398: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
399: const PetscInt eLocal = eLocalx + epr * eLocaly;
401: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
402: }
403: } else if (dim == 3) {
404: const PetscInt epr = stag->nGhost[0];
405: const PetscInt epl = stag->nGhost[0] * stag->nGhost[1];
407: for (PetscInt idx = 0; idx < n; ++idx) {
408: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
409: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
410: const PetscInt eLocalz = pos[idx].k - stag->startGhost[2];
411: const PetscInt eLocal = epl * eLocalz + epr * eLocaly + eLocalx;
413: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
414: }
415: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
416: return 0;
417: }
419: /*@C
420: DMStagVecGetValuesStencil - get vector values using grid indexing
422: Not Collective
424: Input Parameters:
425: + dm - the `DMSTAG` object
426: . vec - the vector object
427: . n - the number of values to obtain
428: - pos - locations to obtain values from (as an array of `DMStagStencil` values)
430: Output Parameter:
431: . val - value at the point
433: Notes:
434: Accepts stencils which refer to global element numbers, but
435: only allows access to entries in the local representation (including ghosts).
437: This approach is not as efficient as getting values directly with `DMStagVecGetArray()`,
438: which is recommended for matrix free operators.
440: Level: advanced
442: .seealso: [](chapter_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecSetValuesStencil()`, `DMStagMatSetValuesStencil()`, `DMStagVecGetArray()`
443: @*/
444: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec, PetscInt n, const DMStagStencil *pos, PetscScalar *val)
445: {
446: DM_Stag *const stag = (DM_Stag *)dm->data;
447: PetscInt nLocal, idx;
448: PetscInt *ix;
449: PetscScalar const *arr;
453: VecGetLocalSize(vec, &nLocal);
455: PetscMalloc1(n, &ix);
456: DMStagStencilToIndexLocal(dm, dm->dim, n, pos, ix);
457: VecGetArrayRead(vec, &arr);
458: for (idx = 0; idx < n; ++idx) val[idx] = arr[ix[idx]];
459: VecRestoreArrayRead(vec, &arr);
460: PetscFree(ix);
461: return 0;
462: }
464: /*@C
465: DMStagVecSetValuesStencil - Set `Vec` values using global grid indexing
467: Not Collective
469: Input Parameters:
470: + dm - the `DMSTAG` object
471: . vec - the `Vec`
472: . n - the number of values to set
473: . pos - the locations to set values, as an array of `DMStagStencil` structs
474: . val - the values to set
475: - insertMode - `INSERT_VALUES` or `ADD_VALUES`
477: Notes:
478: The vector is expected to be a global vector compatible with the DM (usually obtained by `DMGetGlobalVector()` or `DMCreateGlobalVector()`).
480: This approach is not as efficient as setting values directly with `DMStagVecGetArray()`, which is recommended for matrix-free operators.
481: 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()`).
483: Level: advanced
485: .seealso: [](chapter_stag), `DMSTAG`, `Vec`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagMatSetValuesStencil()`, `DMCreateGlobalVector()`, `DMGetLocalVector()`, `DMStagVecGetArray()`
486: @*/
487: PetscErrorCode DMStagVecSetValuesStencil(DM dm, Vec vec, PetscInt n, const DMStagStencil *pos, const PetscScalar *val, InsertMode insertMode)
488: {
489: DM_Stag *const stag = (DM_Stag *)dm->data;
490: PetscInt nLocal;
491: PetscInt *ix;
495: VecGetLocalSize(vec, &nLocal);
497: PetscMalloc1(n, &ix);
498: DMStagStencilToIndexLocal(dm, dm->dim, n, pos, ix);
499: VecSetValuesLocal(vec, n, ix, val, insertMode);
500: PetscFree(ix);
501: return 0;
502: }