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 to 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));
119: PetscFunctionReturn(PETSC_SUCCESS);
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: [](ch_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;
143: PetscFunctionBegin;
145: PetscCall(DMGetDimension(dm, &dim));
146: switch (dim) {
147: case 1:
148: switch (loc) {
149: case DMSTAG_LEFT:
150: case DMSTAG_RIGHT:
151: *dof = stag->dof[0];
152: break;
153: case DMSTAG_ELEMENT:
154: *dof = stag->dof[1];
155: break;
156: default:
157: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
158: }
159: break;
160: case 2:
161: switch (loc) {
162: case DMSTAG_DOWN_LEFT:
163: case DMSTAG_DOWN_RIGHT:
164: case DMSTAG_UP_LEFT:
165: case DMSTAG_UP_RIGHT:
166: *dof = stag->dof[0];
167: break;
168: case DMSTAG_LEFT:
169: case DMSTAG_RIGHT:
170: case DMSTAG_UP:
171: case DMSTAG_DOWN:
172: *dof = stag->dof[1];
173: break;
174: case DMSTAG_ELEMENT:
175: *dof = stag->dof[2];
176: break;
177: default:
178: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
179: }
180: break;
181: case 3:
182: switch (loc) {
183: case DMSTAG_BACK_DOWN_LEFT:
184: case DMSTAG_BACK_DOWN_RIGHT:
185: case DMSTAG_BACK_UP_LEFT:
186: case DMSTAG_BACK_UP_RIGHT:
187: case DMSTAG_FRONT_DOWN_LEFT:
188: case DMSTAG_FRONT_DOWN_RIGHT:
189: case DMSTAG_FRONT_UP_LEFT:
190: case DMSTAG_FRONT_UP_RIGHT:
191: *dof = stag->dof[0];
192: break;
193: case DMSTAG_BACK_DOWN:
194: case DMSTAG_BACK_LEFT:
195: case DMSTAG_BACK_RIGHT:
196: case DMSTAG_BACK_UP:
197: case DMSTAG_DOWN_LEFT:
198: case DMSTAG_DOWN_RIGHT:
199: case DMSTAG_UP_LEFT:
200: case DMSTAG_UP_RIGHT:
201: case DMSTAG_FRONT_DOWN:
202: case DMSTAG_FRONT_LEFT:
203: case DMSTAG_FRONT_RIGHT:
204: case DMSTAG_FRONT_UP:
205: *dof = stag->dof[1];
206: break;
207: case DMSTAG_LEFT:
208: case DMSTAG_RIGHT:
209: case DMSTAG_DOWN:
210: case DMSTAG_UP:
211: case DMSTAG_BACK:
212: case DMSTAG_FRONT:
213: *dof = stag->dof[2];
214: break;
215: case DMSTAG_ELEMENT:
216: *dof = stag->dof[3];
217: break;
218: default:
219: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented for location %s", DMStagStencilLocations[loc]);
220: }
221: break;
222: default:
223: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, dim);
224: }
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*
229: Convert to a location value with only BACK, DOWN, LEFT, and ELEMENT involved
230: */
231: PETSC_INTERN PetscErrorCode DMStagStencilLocationCanonicalize(DMStagStencilLocation loc, DMStagStencilLocation *locCanonical)
232: {
233: PetscFunctionBegin;
234: switch (loc) {
235: case DMSTAG_ELEMENT:
236: *locCanonical = DMSTAG_ELEMENT;
237: break;
238: case DMSTAG_LEFT:
239: case DMSTAG_RIGHT:
240: *locCanonical = DMSTAG_LEFT;
241: break;
242: case DMSTAG_DOWN:
243: case DMSTAG_UP:
244: *locCanonical = DMSTAG_DOWN;
245: break;
246: case DMSTAG_BACK:
247: case DMSTAG_FRONT:
248: *locCanonical = DMSTAG_BACK;
249: break;
250: case DMSTAG_DOWN_LEFT:
251: case DMSTAG_DOWN_RIGHT:
252: case DMSTAG_UP_LEFT:
253: case DMSTAG_UP_RIGHT:
254: *locCanonical = DMSTAG_DOWN_LEFT;
255: break;
256: case DMSTAG_BACK_LEFT:
257: case DMSTAG_BACK_RIGHT:
258: case DMSTAG_FRONT_LEFT:
259: case DMSTAG_FRONT_RIGHT:
260: *locCanonical = DMSTAG_BACK_LEFT;
261: break;
262: case DMSTAG_BACK_DOWN:
263: case DMSTAG_BACK_UP:
264: case DMSTAG_FRONT_DOWN:
265: case DMSTAG_FRONT_UP:
266: *locCanonical = DMSTAG_BACK_DOWN;
267: break;
268: case DMSTAG_BACK_DOWN_LEFT:
269: case DMSTAG_BACK_DOWN_RIGHT:
270: case DMSTAG_BACK_UP_LEFT:
271: case DMSTAG_BACK_UP_RIGHT:
272: case DMSTAG_FRONT_DOWN_LEFT:
273: case DMSTAG_FRONT_DOWN_RIGHT:
274: case DMSTAG_FRONT_UP_LEFT:
275: case DMSTAG_FRONT_UP_RIGHT:
276: *locCanonical = DMSTAG_BACK_DOWN_LEFT;
277: break;
278: default:
279: *locCanonical = DMSTAG_NULL_LOCATION;
280: break;
281: }
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*@C
286: DMStagMatGetValuesStencil - retrieve local matrix entries using grid indexing
288: Not Collective
290: Input Parameters:
291: + dm - the `DMSTAG` object
292: . mat - the matrix
293: . nRow - number of rows
294: . posRow - grid locations (including components) of rows
295: . nCol - number of columns
296: - posCol - grid locations (including components) of columns
298: Output Parameter:
299: . val - logically two-dimensional array of values
301: Level: advanced
303: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagVecSetValuesStencil()`, `DMStagMatSetValuesStencil()`, `MatSetValuesStencil()`, `MatAssemblyBegin()`, `MatAssemblyEnd()`, `DMCreateMatrix()`
304: @*/
305: PetscErrorCode DMStagMatGetValuesStencil(DM dm, Mat mat, PetscInt nRow, const DMStagStencil *posRow, PetscInt nCol, const DMStagStencil *posCol, PetscScalar *val)
306: {
307: PetscInt dim;
308: PetscInt *ir, *ic;
310: PetscFunctionBegin;
313: PetscCall(DMGetDimension(dm, &dim));
314: PetscCall(PetscMalloc2(nRow, &ir, nCol, &ic));
315: PetscCall(DMStagStencilToIndexLocal(dm, dim, nRow, posRow, ir));
316: PetscCall(DMStagStencilToIndexLocal(dm, dim, nCol, posCol, ic));
317: PetscCall(MatGetValuesLocal(mat, nRow, ir, nCol, ic, val));
318: PetscCall(PetscFree2(ir, ic));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /*@C
323: DMStagMatSetValuesStencil - insert or add matrix entries using grid indexing
325: Not Collective
327: Input Parameters:
328: + dm - the `DMSTAG` object
329: . mat - the matrix
330: . nRow - number of rows
331: . posRow - grid locations (including components) of rows
332: . nCol - number of columns
333: . posCol - grid locations (including components) of columns
334: . val - logically two-dimensional array of values
335: - insertMode - `INSERT_VALUES` or `ADD_VALUES`
337: Notes:
338: See notes for `MatSetValuesStencil()`
340: Level: intermediate
342: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagVecSetValuesStencil()`, `DMStagMatGetValuesStencil()`, `MatSetValuesStencil()`, `MatAssemblyBegin()`, `MatAssemblyEnd()`, `DMCreateMatrix()`
343: @*/
344: PetscErrorCode DMStagMatSetValuesStencil(DM dm, Mat mat, PetscInt nRow, const DMStagStencil *posRow, PetscInt nCol, const DMStagStencil *posCol, const PetscScalar *val, InsertMode insertMode)
345: {
346: PetscInt *ir, *ic;
348: PetscFunctionBegin;
351: PetscCall(PetscMalloc2(nRow, &ir, nCol, &ic));
352: PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, nRow, posRow, ir));
353: PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, nCol, posCol, ic));
354: PetscCall(MatSetValuesLocal(mat, nRow, ir, nCol, ic, val, insertMode));
355: PetscCall(PetscFree2(ir, ic));
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: /*@C
360: DMStagStencilToIndexLocal - Convert an array of `DMStagStenci`l objects to an array of indices into a local vector.
362: Not Collective
364: Input Parameters:
365: + dm - the `DMSTAG` object
366: . dim - the dimension of the `DMSTAG` object
367: . n - the number of `DMStagStencil` objects
368: - pos - an array of `n` `DMStagStencil` objects
370: Output Parameter:
371: . ix - output array of `n` indices
373: Notes:
374: The `DMStagStencil` objects in `pos` use global element indices.
376: The `.c` fields in `pos` must always be set (even if to `0`).
378: Developer Notes:
379: This is a "hot" function, and accepts the dimension redundantly to avoid having to perform any error checking inside the function.
381: Level: developer
383: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencilLocation`, `DMStagStencil`, `DMGetLocalVector`, `DMCreateLocalVector`
384: @*/
385: PetscErrorCode DMStagStencilToIndexLocal(DM dm, PetscInt dim, PetscInt n, const DMStagStencil *pos, PetscInt *ix)
386: {
387: const DM_Stag *const stag = (DM_Stag *)dm->data;
388: const PetscInt epe = stag->entriesPerElement;
390: PetscFunctionBeginHot;
391: if (dim == 1) {
392: for (PetscInt idx = 0; idx < n; ++idx) {
393: const PetscInt eLocal = pos[idx].i - stag->startGhost[0];
395: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
396: }
397: } else if (dim == 2) {
398: const PetscInt epr = stag->nGhost[0];
400: for (PetscInt idx = 0; idx < n; ++idx) {
401: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
402: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
403: const PetscInt eLocal = eLocalx + epr * eLocaly;
405: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
406: }
407: } else if (dim == 3) {
408: const PetscInt epr = stag->nGhost[0];
409: const PetscInt epl = stag->nGhost[0] * stag->nGhost[1];
411: for (PetscInt idx = 0; idx < n; ++idx) {
412: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
413: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
414: const PetscInt eLocalz = pos[idx].k - stag->startGhost[2];
415: const PetscInt eLocal = epl * eLocalz + epr * eLocaly + eLocalx;
417: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
418: }
419: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %" PetscInt_FMT, dim);
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: /*@C
424: DMStagVecGetValuesStencil - get vector values using grid indexing
426: Not Collective
428: Input Parameters:
429: + dm - the `DMSTAG` object
430: . vec - the vector object
431: . n - the number of values to obtain
432: - pos - locations to obtain values from (as an array of `DMStagStencil` values)
434: Output Parameter:
435: . val - value at the point
437: Notes:
438: Accepts stencils which refer to global element numbers, but
439: only allows access to entries in the local representation (including ghosts).
441: This approach is not as efficient as getting values directly with `DMStagVecGetArray()`,
442: which is recommended for matrix-free operators.
444: Level: advanced
446: .seealso: [](ch_stag), `DMSTAG`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecSetValuesStencil()`, `DMStagMatSetValuesStencil()`, `DMStagVecGetArray()`
447: @*/
448: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec, PetscInt n, const DMStagStencil *pos, PetscScalar *val)
449: {
450: DM_Stag *const stag = (DM_Stag *)dm->data;
451: PetscInt nLocal, idx;
452: PetscInt *ix;
453: PetscScalar const *arr;
455: PetscFunctionBegin;
458: PetscCall(VecGetLocalSize(vec, &nLocal));
459: 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);
460: PetscCall(PetscMalloc1(n, &ix));
461: PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, n, pos, ix));
462: PetscCall(VecGetArrayRead(vec, &arr));
463: for (idx = 0; idx < n; ++idx) val[idx] = arr[ix[idx]];
464: PetscCall(VecRestoreArrayRead(vec, &arr));
465: PetscCall(PetscFree(ix));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@C
470: DMStagVecSetValuesStencil - Set `Vec` values using global grid indexing
472: Not Collective
474: Input Parameters:
475: + dm - the `DMSTAG` object
476: . vec - the `Vec`
477: . n - the number of values to set
478: . pos - the locations to set values, as an array of `DMStagStencil` structs
479: . val - the values to set
480: - insertMode - `INSERT_VALUES` or `ADD_VALUES`
482: Notes:
483: The vector is expected to be a global vector compatible with the DM (usually obtained by `DMGetGlobalVector()` or `DMCreateGlobalVector()`).
485: This approach is not as efficient as setting values directly with `DMStagVecGetArray()`, which is recommended for matrix-free operators.
486: 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()`).
488: Level: advanced
490: .seealso: [](ch_stag), `DMSTAG`, `Vec`, `DMStagStencil`, `DMStagStencilLocation`, `DMStagVecGetValuesStencil()`, `DMStagMatSetValuesStencil()`, `DMCreateGlobalVector()`, `DMGetLocalVector()`, `DMStagVecGetArray()`
491: @*/
492: PetscErrorCode DMStagVecSetValuesStencil(DM dm, Vec vec, PetscInt n, const DMStagStencil *pos, const PetscScalar *val, InsertMode insertMode)
493: {
494: DM_Stag *const stag = (DM_Stag *)dm->data;
495: PetscInt nLocal;
496: PetscInt *ix;
498: PetscFunctionBegin;
501: PetscCall(VecGetLocalSize(vec, &nLocal));
502: 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);
503: PetscCall(PetscMalloc1(n, &ix));
504: PetscCall(DMStagStencilToIndexLocal(dm, dm->dim, n, pos, ix));
505: PetscCall(VecSetValuesLocal(vec, n, ix, val, insertMode));
506: PetscCall(PetscFree(ix));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }