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, &ltog);
 79:   ISLocalToGlobalMappingGetIndices(ltog, &ltogidx);
 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, &ltogidx);
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: }