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: [](ch_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: PetscFunctionBegin;
41: PetscCall(DMGetDimension(dm, &dim));
42: PetscCheck(dim >= 1 && dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
44: /* To assert that the resulting IS has unique, sorted, entries, we perform
45: a bucket sort, taking advantage of the fact that DMStagStencilLocation
46: enum values are integers starting with 1, in canonical order */
47: nc_max = 1; // maximum number of components to represent these stencils
48: n_stencil_unique = 0;
49: for (PetscInt p = 0; p < n_stencil; ++p) nc_max = PetscMax(nc_max, (stencils[p].c + 1));
50: PetscCall(PetscCalloc1(DMSTAG_NUMBER_LOCATIONS * nc_max, &stencil_active));
51: for (PetscInt p = 0; p < n_stencil; ++p) {
52: DMStagStencilLocation loc_canonical;
53: PetscInt slot;
55: PetscCall(DMStagStencilLocationCanonicalize(stencils[p].loc, &loc_canonical));
56: slot = nc_max * ((PetscInt)loc_canonical) + stencils[p].c;
57: if (stencil_active[slot] == 0) {
58: stencil_active[slot] = 1;
59: ++n_stencil_unique;
60: }
61: }
62: PetscCall(PetscMalloc1(n_stencil_unique, &stencils_ordered_unique));
63: {
64: PetscInt p = 0;
66: for (PetscInt i = 1; i < DMSTAG_NUMBER_LOCATIONS; ++i) {
67: for (PetscInt c = 0; c < nc_max; ++c) {
68: if (stencil_active[nc_max * i + c] != 0) {
69: stencils_ordered_unique[p].loc = (DMStagStencilLocation)i;
70: stencils_ordered_unique[p].c = c;
71: ++p;
72: }
73: }
74: }
75: }
76: PetscCall(PetscFree(stencil_active));
78: PetscCall(PetscMalloc1(n_stencil_unique, &idxLocal));
79: PetscCall(DMGetLocalToGlobalMapping(dm, <og));
80: PetscCall(ISLocalToGlobalMappingGetIndices(ltog, <ogidx));
81: PetscCall(DMStagGetCorners(dm, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &extraPoint[0], &extraPoint[1], &extraPoint[2]));
82: for (PetscInt d = dim; d < DMSTAG_MAX_DIM; ++d) {
83: start[d] = 0;
84: n[d] = 1; /* To allow for a single loop nest below */
85: extraPoint[d] = 0;
86: }
87: nidx = n_stencil_unique;
88: for (PetscInt d = 0; d < dim; ++d) nidx *= (n[d] + 1); /* Overestimate (always assumes extraPoint) */
89: PetscCall(PetscMalloc1(nidx, &idx));
90: count = 0;
91: /* Note that unused loop variables are not accessed, for lower dimensions */
92: for (PetscInt k = start[2]; k < start[2] + n[2] + extraPoint[2]; ++k) {
93: for (PetscInt j = start[1]; j < start[1] + n[1] + extraPoint[1]; ++j) {
94: for (PetscInt i = start[0]; i < start[0] + n[0] + extraPoint[0]; ++i) {
95: for (PetscInt p = 0; p < n_stencil_unique; ++p) {
96: stencils_ordered_unique[p].i = i;
97: stencils_ordered_unique[p].j = j;
98: stencils_ordered_unique[p].k = k;
99: }
100: PetscCall(DMStagStencilToIndexLocal(dm, dim, n_stencil_unique, stencils_ordered_unique, idxLocal));
101: for (PetscInt p = 0; p < n_stencil_unique; ++p) {
102: const PetscInt gidx = ltogidx[idxLocal[p]];
103: if (gidx >= 0) {
104: idx[count] = gidx;
105: ++count;
106: }
107: }
108: }
109: }
110: }
112: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltog, <ogidx));
113: PetscCall(PetscFree(stencils_ordered_unique));
114: PetscCall(PetscFree(idxLocal));
116: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), count, idx, PETSC_OWN_POINTER, is));
117: PetscCall(ISSetInfo(*is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
118: PetscCall(ISSetInfo(*is, IS_UNIQUE, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: /*@C
124: DMStagGetLocationDOF - Get number of DOF associated with a given point in a `DMSTAG` grid
126: Not Collective
128: Input Parameters:
129: + dm - the `DMSTAG` object
130: - loc - grid point (see `DMStagStencilLocation`)
132: Output Parameter:
133: . dof - the number of DOF (components) living at `loc` in `dm`
135: Level: intermediate
137: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencilLocation`, `DMStagStencil`, `DMDAGetDof()`
138: @*/
139: PetscErrorCode DMStagGetLocationDOF(DM dm, DMStagStencilLocation loc, PetscInt *dof)
140: {
141: const DM_Stag *const stag = (DM_Stag *)dm->data;
142: PetscInt dim;
144: PetscFunctionBegin;
146: PetscCall(DMGetDimension(dm, &dim));
147: switch (dim) {
148: case 1:
149: switch (loc) {
150: case DMSTAG_LEFT:
151: case DMSTAG_RIGHT:
152: *dof = stag->dof[0];
153: break;
154: case DMSTAG_ELEMENT:
155: *dof = stag->dof[1];
156: break;
157: default:
158: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
159: }
160: break;
161: case 2:
162: switch (loc) {
163: case DMSTAG_DOWN_LEFT:
164: case DMSTAG_DOWN_RIGHT:
165: case DMSTAG_UP_LEFT:
166: case DMSTAG_UP_RIGHT:
167: *dof = stag->dof[0];
168: break;
169: case DMSTAG_LEFT:
170: case DMSTAG_RIGHT:
171: case DMSTAG_UP:
172: case DMSTAG_DOWN:
173: *dof = stag->dof[1];
174: break;
175: case DMSTAG_ELEMENT:
176: *dof = stag->dof[2];
177: break;
178: default:
179: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
180: }
181: break;
182: case 3:
183: switch (loc) {
184: case DMSTAG_BACK_DOWN_LEFT:
185: case DMSTAG_BACK_DOWN_RIGHT:
186: case DMSTAG_BACK_UP_LEFT:
187: case DMSTAG_BACK_UP_RIGHT:
188: case DMSTAG_FRONT_DOWN_LEFT:
189: case DMSTAG_FRONT_DOWN_RIGHT:
190: case DMSTAG_FRONT_UP_LEFT:
191: case DMSTAG_FRONT_UP_RIGHT:
192: *dof = stag->dof[0];
193: break;
194: case DMSTAG_BACK_DOWN:
195: case DMSTAG_BACK_LEFT:
196: case DMSTAG_BACK_RIGHT:
197: case DMSTAG_BACK_UP:
198: case DMSTAG_DOWN_LEFT:
199: case DMSTAG_DOWN_RIGHT:
200: case DMSTAG_UP_LEFT:
201: case DMSTAG_UP_RIGHT:
202: case DMSTAG_FRONT_DOWN:
203: case DMSTAG_FRONT_LEFT:
204: case DMSTAG_FRONT_RIGHT:
205: case DMSTAG_FRONT_UP:
206: *dof = stag->dof[1];
207: break;
208: case DMSTAG_LEFT:
209: case DMSTAG_RIGHT:
210: case DMSTAG_DOWN:
211: case DMSTAG_UP:
212: case DMSTAG_BACK:
213: case DMSTAG_FRONT:
214: *dof = stag->dof[2];
215: break;
216: case DMSTAG_ELEMENT:
217: *dof = stag->dof[3];
218: break;
219: default:
220: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
221: }
222: break;
223: default:
224: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
225: }
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*
230: Convert to a location value with only BACK, DOWN, LEFT, and ELEMENT involved
231: */
232: PETSC_INTERN PetscErrorCode DMStagStencilLocationCanonicalize(DMStagStencilLocation loc, DMStagStencilLocation *locCanonical)
233: {
234: PetscFunctionBegin;
235: switch (loc) {
236: case DMSTAG_ELEMENT:
237: *locCanonical = DMSTAG_ELEMENT;
238: break;
239: case DMSTAG_LEFT:
240: case DMSTAG_RIGHT:
241: *locCanonical = DMSTAG_LEFT;
242: break;
243: case DMSTAG_DOWN:
244: case DMSTAG_UP:
245: *locCanonical = DMSTAG_DOWN;
246: break;
247: case DMSTAG_BACK:
248: case DMSTAG_FRONT:
249: *locCanonical = DMSTAG_BACK;
250: break;
251: case DMSTAG_DOWN_LEFT:
252: case DMSTAG_DOWN_RIGHT:
253: case DMSTAG_UP_LEFT:
254: case DMSTAG_UP_RIGHT:
255: *locCanonical = DMSTAG_DOWN_LEFT;
256: break;
257: case DMSTAG_BACK_LEFT:
258: case DMSTAG_BACK_RIGHT:
259: case DMSTAG_FRONT_LEFT:
260: case DMSTAG_FRONT_RIGHT:
261: *locCanonical = DMSTAG_BACK_LEFT;
262: break;
263: case DMSTAG_BACK_DOWN:
264: case DMSTAG_BACK_UP:
265: case DMSTAG_FRONT_DOWN:
266: case DMSTAG_FRONT_UP:
267: *locCanonical = DMSTAG_BACK_DOWN;
268: break;
269: case DMSTAG_BACK_DOWN_LEFT:
270: case DMSTAG_BACK_DOWN_RIGHT:
271: case DMSTAG_BACK_UP_LEFT:
272: case DMSTAG_BACK_UP_RIGHT:
273: case DMSTAG_FRONT_DOWN_LEFT:
274: case DMSTAG_FRONT_DOWN_RIGHT:
275: case DMSTAG_FRONT_UP_LEFT:
276: case DMSTAG_FRONT_UP_RIGHT:
277: *locCanonical = DMSTAG_BACK_DOWN_LEFT;
278: break;
279: default:
280: *locCanonical = DMSTAG_NULL_LOCATION;
281: break;
282: }
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@C
287: DMStagMatGetValuesStencil - retrieve local matrix entries using grid indexing
289: Not Collective
291: Input Parameters:
292: + dm - the `DMSTAG` object
293: . mat - the matrix
294: . nRow - number of rows
295: . posRow - grid locations (including components) of rows
296: . nCol - number of columns
297: - posCol - grid locations (including components) of columns
299: Output Parameter:
300: . val - logically two-dimensional array of values
302: Level: advanced
304: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagVecSetValuesStencil()`, `DMStagMatSetValuesStencil()`, `MatSetValuesStencil()`, `MatAssemblyBegin()`, `MatAssemblyEnd()`, `DMCreateMatrix()`
305: @*/
306: PetscErrorCode DMStagMatGetValuesStencil(DM dm, Mat mat, PetscInt nRow, const DMStagStencil *posRow, PetscInt nCol, const DMStagStencil *posCol, PetscScalar *val)
307: {
308: PetscInt dim;
309: PetscInt *ir, *ic;
311: PetscFunctionBegin;
314: PetscCall(DMGetDimension(dm, &dim));
315: PetscCall(PetscMalloc2(nRow, &ir, nCol, &ic));
316: PetscCall(DMStagStencilToIndexLocal(dm, dim, nRow, posRow, ir));
317: PetscCall(DMStagStencilToIndexLocal(dm, dim, nCol, posCol, ic));
318: PetscCall(MatGetValuesLocal(mat, nRow, ir, nCol, ic, val));
319: PetscCall(PetscFree2(ir, ic));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /*@C
324: DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing
326: Not Collective
328: Input Parameters:
329: + dm - the `DMSTAG` object
330: . mat - the matrix
331: . nRow - number of rows
332: . posRow - grid locations (including components) of rows
333: . nCol - number of columns
334: . posCol - grid locations (including components) of columns
335: . val - logically two-dimensional array of values
336: - insertMode - `INSERT_VALUES` or `ADD_VALUES`
338: Notes:
339: See notes for `MatSetValuesStencil()`
341: Level: intermediate
343: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagVecSetValuesStencil()`, `DMStagMatGetValuesStencil()`, `MatSetValuesStencil()`, `MatAssemblyBegin()`, `MatAssemblyEnd()`, `DMCreateMatrix()`
344: @*/
345: PetscErrorCode DMStagMatSetValuesStencil(DM dm, Mat mat, PetscInt nRow, const DMStagStencil *posRow, PetscInt nCol, const DMStagStencil *posCol, const PetscScalar *val, InsertMode insertMode)
346: {
347: PetscInt *ir, *ic;
349: PetscFunctionBegin;
352: PetscCall(PetscMalloc2(nRow, &ir, nCol, &ic));
353: PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, nRow, posRow, ir));
354: PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, nCol, posCol, ic));
355: PetscCall(MatSetValuesLocal(mat, nRow, ir, nCol, ic, val, insertMode));
356: PetscCall(PetscFree2(ir, ic));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: /*@C
361: DMStagStencilToIndexLocal - Convert an array of `DMStagStenci`l objects to an array of indices into a local vector.
363: Not Collective
365: Input Parameters:
366: + dm - the `DMSTAG` object
367: . dim - the dimension of the `DMSTAG` object
368: . n - the number of `DMStagStencil` objects
369: - pos - an array of `n` `DMStagStencil` objects
371: Output Parameter:
372: . ix - output array of `n` indices
374: Notes:
375: The `DMStagStencil` objects in `pos` use global element indices.
377: The `.c` fields in `pos` must always be set (even if to `0`).
379: Developer Notes:
380: This is a "hot" function, and accepts the dimension redundantly to avoid having to perform any error checking inside the function.
382: Level: developer
384: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencilLocation`, `DMStagStencil`, `DMGetLocalVector`, `DMCreateLocalVector`
385: @*/
386: PetscErrorCode DMStagStencilToIndexLocal(DM dm, PetscInt dim, PetscInt n, const DMStagStencil *pos, PetscInt *ix)
387: {
388: const DM_Stag *const stag = (DM_Stag *)dm->data;
389: const PetscInt epe = stag->entriesPerElement;
391: PetscFunctionBeginHot;
392: if (dim == 1) {
393: for (PetscInt idx = 0; idx < n; ++idx) {
394: const PetscInt eLocal = pos[idx].i - stag->startGhost[0];
396: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
397: }
398: } else if (dim == 2) {
399: const PetscInt epr = stag->nGhost[0];
401: for (PetscInt idx = 0; idx < n; ++idx) {
402: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
403: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
404: const PetscInt eLocal = eLocalx + epr * eLocaly;
406: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
407: }
408: } else if (dim == 3) {
409: const PetscInt epr = stag->nGhost[0];
410: const PetscInt epl = stag->nGhost[0] * stag->nGhost[1];
412: for (PetscInt idx = 0; idx < n; ++idx) {
413: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
414: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
415: const PetscInt eLocalz = pos[idx].k - stag->startGhost[2];
416: const PetscInt eLocal = epl * eLocalz + epr * eLocaly + eLocalx;
418: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
419: }
420: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@C
425: DMStagVecGetValuesStencil - get vector values using grid indexing
427: Not Collective
429: Input Parameters:
430: + dm - the `DMSTAG` object
431: . vec - the vector object
432: . n - the number of values to obtain
433: - pos - locations to obtain values from (as an array of `DMStagStencil` values)
435: Output Parameter:
436: . val - value at the point
438: Notes:
439: Accepts stencils which refer to global element numbers, but
440: only allows access to entries in the local representation (including ghosts).
442: This approach is not as efficient as getting values directly with `DMStagVecGetArray()`,
443: which is recommended for matrix-free operators.
445: Level: advanced
447: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecSetValuesStencil()`, `DMStagMatSetValuesStencil()`, `DMStagVecGetArray()`
448: @*/
449: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec, PetscInt n, const DMStagStencil *pos, PetscScalar *val)
450: {
451: DM_Stag *const stag = (DM_Stag *)dm->data;
452: PetscInt nLocal, idx;
453: PetscInt *ix;
454: PetscScalar const *arr;
456: PetscFunctionBegin;
459: PetscCall(VecGetLocalSize(vec, &nLocal));
460: PetscCheck(nLocal == stag->entriesGhost, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector should be a local vector. Local size %" PetscInt_FMT " does not match expected %" PetscInt_FMT, nLocal, stag->entriesGhost);
461: PetscCall(PetscMalloc1(n, &ix));
462: PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, n, pos, ix));
463: PetscCall(VecGetArrayRead(vec, &arr));
464: for (idx = 0; idx < n; ++idx) val[idx] = arr[ix[idx]];
465: PetscCall(VecRestoreArrayRead(vec, &arr));
466: PetscCall(PetscFree(ix));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*@C
471: DMStagVecSetValuesStencil - Set `Vec` values using global grid indexing
473: Not Collective
475: Input Parameters:
476: + dm - the `DMSTAG` object
477: . vec - the `Vec`
478: . n - the number of values to set
479: . pos - the locations to set values, as an array of `DMStagStencil` structs
480: . val - the values to set
481: - insertMode - `INSERT_VALUES` or `ADD_VALUES`
483: Notes:
484: The vector is expected to be a global vector compatible with the DM (usually obtained by `DMGetGlobalVector()` or `DMCreateGlobalVector()`).
486: This approach is not as efficient as setting values directly with `DMStagVecGetArray()`, which is recommended for matrix-free operators.
487: 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()`).
489: Level: advanced
491: .seealso: [](ch_stag), `DMSTAG`, `Vec`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagMatSetValuesStencil()`, `DMCreateGlobalVector()`, `DMGetLocalVector()`, `DMStagVecGetArray()`
492: @*/
493: PetscErrorCode DMStagVecSetValuesStencil(DM dm, Vec vec, PetscInt n, const DMStagStencil *pos, const PetscScalar *val, InsertMode insertMode)
494: {
495: DM_Stag *const stag = (DM_Stag *)dm->data;
496: PetscInt nLocal;
497: PetscInt *ix;
499: PetscFunctionBegin;
502: PetscCall(VecGetLocalSize(vec, &nLocal));
503: PetscCheck(nLocal == stag->entries, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Provided vec has a different number of local entries (%" PetscInt_FMT ") than expected (%" PetscInt_FMT "). It should be a global vector", nLocal, stag->entries);
504: PetscCall(PetscMalloc1(n, &ix));
505: PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, n, pos, ix));
506: PetscCall(VecSetValuesLocal(vec, n, ix, val, insertMode));
507: PetscCall(PetscFree(ix));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }