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: . nStencil - 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 s are ignored
26: Level: advanced
28: .seealso: DMSTAG, IS, DMStagStencil, DMCreateGlobalVector
29: @*/
30: PetscErrorCode DMStagCreateISFromStencils(DM dm,PetscInt nStencil,DMStagStencil* stencils,IS *is)
31: {
32: DMStagStencil *ss;
33: PetscInt *idx,*idxLocal;
34: const PetscInt *ltogidx;
35: PetscInt p,p2,pmax,i,j,k,d,dim,count,nidx;
36: ISLocalToGlobalMapping ltog;
37: PetscInt start[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],extraPoint[DMSTAG_MAX_DIM];
39: DMGetDimension(dm,&dim);
42: /* Only use non-redundant stencils */
43: PetscMalloc1(nStencil,&ss);
44: pmax = 0;
45: for (p=0; p<nStencil; ++p) {
46: PetscBool skip = PETSC_FALSE;
47: DMStagStencil stencilPotential = stencils[p];
48: DMStagStencilLocationCanonicalize(stencils[p].loc,&stencilPotential.loc);
49: for (p2=0; p2<pmax; ++p2) { /* Quadratic complexity algorithm in nStencil */
50: if (stencilPotential.loc == ss[p2].loc && stencilPotential.c == ss[p2].c) {
51: skip = PETSC_TRUE;
52: break;
53: }
54: }
55: if (!skip) {
56: ss[pmax] = stencilPotential;
57: ++pmax;
58: }
59: }
61: PetscMalloc1(pmax,&idxLocal);
62: DMGetLocalToGlobalMapping(dm,<og);
63: ISLocalToGlobalMappingGetIndices(ltog,<ogidx);
64: DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&extraPoint[0],&extraPoint[1],&extraPoint[2]);
65: for (d=dim; d<DMSTAG_MAX_DIM; ++d) {
66: start[d] = 0;
67: n[d] = 1; /* To allow for a single loop nest below */
68: extraPoint[d] = 0;
69: }
70: nidx = pmax; for (d=0; d<dim; ++d) nidx *= (n[d]+1); /* Overestimate (always assumes extraPoint) */
71: PetscMalloc1(nidx,&idx);
72: count = 0;
73: /* Note that unused loop variables are not accessed, for lower dimensions */
74: for (k=start[2]; k<start[2]+n[2]+extraPoint[2]; ++k) {
75: for (j=start[1]; j<start[1]+n[1]+extraPoint[1]; ++j) {
76: for (i=start[0]; i<start[0]+n[0]+extraPoint[0]; ++i) {
77: for (p=0; p<pmax; ++p) {
78: ss[p].i = i; ss[p].j = j; ss[p].k = k;
79: }
80: DMStagStencilToIndexLocal(dm,dim,pmax,ss,idxLocal);
81: for (p=0; p<pmax; ++p) {
82: const PetscInt gidx = ltogidx[idxLocal[p]];
83: if (gidx >= 0) {
84: idx[count] = gidx;
85: ++count;
86: }
87: }
88: }
89: }
90: }
91: ISLocalToGlobalMappingRestoreIndices(ltog,<ogidx);
92: ISCreateGeneral(PetscObjectComm((PetscObject)dm),count,idx,PETSC_OWN_POINTER,is);
94: PetscFree(ss);
95: PetscFree(idxLocal);
96: return 0;
97: }
99: /*@C
100: DMStagGetLocationDOF - Get number of DOF associated with a given point in a DMStag grid
102: Not Collective
104: Input Parameters:
105: + dm - the DMStag object
106: - loc - grid point (see DMStagStencilLocation)
108: Output Parameter:
109: . dof - the number of dof (components) living at loc in dm
111: Level: intermediate
113: .seealso: DMSTAG, DMStagStencilLocation, DMStagStencil, DMDAGetDof()
114: @*/
115: PetscErrorCode DMStagGetLocationDOF(DM dm,DMStagStencilLocation loc,PetscInt *dof)
116: {
117: const DM_Stag * const stag = (DM_Stag*)dm->data;
118: PetscInt dim;
121: DMGetDimension(dm,&dim);
122: switch (dim) {
123: case 1:
124: switch (loc) {
125: case DMSTAG_LEFT:
126: case DMSTAG_RIGHT:
127: *dof = stag->dof[0]; break;
128: case DMSTAG_ELEMENT:
129: *dof = stag->dof[1]; break;
130: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location %s",DMStagStencilLocations[loc]);
131: }
132: break;
133: case 2:
134: switch (loc) {
135: case DMSTAG_DOWN_LEFT:
136: case DMSTAG_DOWN_RIGHT:
137: case DMSTAG_UP_LEFT:
138: case DMSTAG_UP_RIGHT:
139: *dof = stag->dof[0]; break;
140: case DMSTAG_LEFT:
141: case DMSTAG_RIGHT:
142: case DMSTAG_UP:
143: case DMSTAG_DOWN:
144: *dof = stag->dof[1]; break;
145: case DMSTAG_ELEMENT:
146: *dof = stag->dof[2]; break;
147: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location %s",DMStagStencilLocations[loc]);
148: }
149: break;
150: case 3:
151: switch (loc) {
152: case DMSTAG_BACK_DOWN_LEFT:
153: case DMSTAG_BACK_DOWN_RIGHT:
154: case DMSTAG_BACK_UP_LEFT:
155: case DMSTAG_BACK_UP_RIGHT:
156: case DMSTAG_FRONT_DOWN_LEFT:
157: case DMSTAG_FRONT_DOWN_RIGHT:
158: case DMSTAG_FRONT_UP_LEFT:
159: case DMSTAG_FRONT_UP_RIGHT:
160: *dof = stag->dof[0]; break;
161: case DMSTAG_BACK_DOWN:
162: case DMSTAG_BACK_LEFT:
163: case DMSTAG_BACK_RIGHT:
164: case DMSTAG_BACK_UP:
165: case DMSTAG_DOWN_LEFT:
166: case DMSTAG_DOWN_RIGHT:
167: case DMSTAG_UP_LEFT:
168: case DMSTAG_UP_RIGHT:
169: case DMSTAG_FRONT_DOWN:
170: case DMSTAG_FRONT_LEFT:
171: case DMSTAG_FRONT_RIGHT:
172: case DMSTAG_FRONT_UP:
173: *dof = stag->dof[1]; break;
174: case DMSTAG_LEFT:
175: case DMSTAG_RIGHT:
176: case DMSTAG_DOWN:
177: case DMSTAG_UP:
178: case DMSTAG_BACK:
179: case DMSTAG_FRONT:
180: *dof = stag->dof[2]; break;
181: case DMSTAG_ELEMENT:
182: *dof = stag->dof[3]; break;
183: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for location %s",DMStagStencilLocations[loc]);
184: }
185: break;
186: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dimension %D",dim);
187: }
188: return 0;
189: }
191: /*
192: Convert to a location value with only BACK, DOWN, LEFT, and ELEMENT involved
193: */
194: PETSC_INTERN PetscErrorCode DMStagStencilLocationCanonicalize(DMStagStencilLocation loc,DMStagStencilLocation *locCanonical)
195: {
196: switch (loc) {
197: case DMSTAG_ELEMENT:
198: *locCanonical = DMSTAG_ELEMENT;
199: break;
200: case DMSTAG_LEFT:
201: case DMSTAG_RIGHT:
202: *locCanonical = DMSTAG_LEFT;
203: break;
204: case DMSTAG_DOWN:
205: case DMSTAG_UP:
206: *locCanonical = DMSTAG_DOWN;
207: break;
208: case DMSTAG_BACK:
209: case DMSTAG_FRONT:
210: *locCanonical = DMSTAG_BACK;
211: break;
212: case DMSTAG_DOWN_LEFT :
213: case DMSTAG_DOWN_RIGHT :
214: case DMSTAG_UP_LEFT :
215: case DMSTAG_UP_RIGHT :
216: *locCanonical = DMSTAG_DOWN_LEFT;
217: break;
218: case DMSTAG_BACK_LEFT:
219: case DMSTAG_BACK_RIGHT:
220: case DMSTAG_FRONT_LEFT:
221: case DMSTAG_FRONT_RIGHT:
222: *locCanonical = DMSTAG_BACK_LEFT;
223: break;
224: case DMSTAG_BACK_DOWN:
225: case DMSTAG_BACK_UP:
226: case DMSTAG_FRONT_DOWN:
227: case DMSTAG_FRONT_UP:
228: *locCanonical = DMSTAG_BACK_DOWN;
229: break;
230: case DMSTAG_BACK_DOWN_LEFT:
231: case DMSTAG_BACK_DOWN_RIGHT:
232: case DMSTAG_BACK_UP_LEFT:
233: case DMSTAG_BACK_UP_RIGHT:
234: case DMSTAG_FRONT_DOWN_LEFT:
235: case DMSTAG_FRONT_DOWN_RIGHT:
236: case DMSTAG_FRONT_UP_LEFT:
237: case DMSTAG_FRONT_UP_RIGHT:
238: *locCanonical = DMSTAG_BACK_DOWN_LEFT;
239: break;
240: default :
241: *locCanonical = DMSTAG_NULL_LOCATION;
242: break;
243: }
244: return 0;
245: }
247: /*@C
248: DMStagMatGetValuesStencil - retrieve local matrix entries using grid indexing
250: Not Collective
252: Input Parameters:
253: + dm - the DMStag object
254: . mat - the matrix
255: . nRow - number of rows
256: . posRow - grid locations (including components) of rows
257: . nCol - number of columns
258: - posCol - grid locations (including components) of columns
260: Output Parameter:
261: . val - logically two-dimensional array of values
263: Level: advanced
265: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
266: @*/
267: PetscErrorCode DMStagMatGetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,PetscScalar *val)
268: {
269: PetscInt dim;
270: PetscInt *ir,*ic;
274: DMGetDimension(dm,&dim);
275: PetscMalloc2(nRow,&ir,nCol,&ic);
276: DMStagStencilToIndexLocal(dm,dim,nRow,posRow,ir);
277: DMStagStencilToIndexLocal(dm,dim,nCol,posCol,ic);
278: MatGetValuesLocal(mat,nRow,ir,nCol,ic,val);
279: PetscFree2(ir,ic);
280: return 0;
281: }
283: /*@C
284: DMStagMatSetValuesStencil - insert or add 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
295: . val - logically two-dimensional array of values
296: - insertmode - INSERT_VALUES or ADD_VALUES
298: Notes:
299: See notes for MatSetValuesStencil()
301: Level: intermediate
303: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagVecSetValuesStencil(), DMStagMatGetValuesStencil(), MatSetValuesStencil(), MatAssemblyBegin(), MatAssemblyEnd(), DMCreateMatrix()
304: @*/
305: PetscErrorCode DMStagMatSetValuesStencil(DM dm,Mat mat,PetscInt nRow,const DMStagStencil *posRow,PetscInt nCol,const DMStagStencil *posCol,const PetscScalar *val,InsertMode insertMode)
306: {
307: PetscInt *ir,*ic;
311: PetscMalloc2(nRow,&ir,nCol,&ic);
312: DMStagStencilToIndexLocal(dm,dm->dim,nRow,posRow,ir);
313: DMStagStencilToIndexLocal(dm,dm->dim,nCol,posCol,ic);
314: MatSetValuesLocal(mat,nRow,ir,nCol,ic,val,insertMode);
315: PetscFree2(ir,ic);
316: return 0;
317: }
319: /*@C
320: DMStagStencilToIndexLocal - Convert an array of DMStagStencil objects to an array of indices into a local vector.
322: Not Collective
324: Input Parameters:
325: + dm - the DMStag object
326: . dim - the dimension of the DMStag object
327: . n - the number of DMStagStencil objects
328: - pos - an array of n DMStagStencil objects
330: Output Parameter:
331: . ix - output array of n indices
333: Notes:
334: The DMStagStencil objects in pos use global element indices.
336: The .c fields in pos must always be set (even if to 0).
338: Developer Notes:
339: This is a "hot" function, and accepts the dimension redundantly to avoid having to perform any error checking inside the function.
341: Level: developer
343: .seealso: DMSTAG, DMStagStencilLocation, DMStagStencil, DMGetLocalVector, DMCreateLocalVector
344: @*/
345: PetscErrorCode DMStagStencilToIndexLocal(DM dm,PetscInt dim,PetscInt n,const DMStagStencil *pos,PetscInt *ix)
346: {
347: const DM_Stag * const stag = (DM_Stag*)dm->data;
348: const PetscInt epe = stag->entriesPerElement;
351: if (dim == 1) {
352: for (PetscInt idx=0; idx<n; ++idx) {
353: const PetscInt eLocal = pos[idx].i - stag->startGhost[0];
355: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
356: }
357: } else if (dim == 2) {
358: const PetscInt epr = stag->nGhost[0];
360: for (PetscInt idx=0; idx<n; ++idx) {
361: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
362: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
363: const PetscInt eLocal = eLocalx + epr*eLocaly;
365: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
366: }
367: } else if (dim == 3) {
368: const PetscInt epr = stag->nGhost[0];
369: const PetscInt epl = stag->nGhost[0]*stag->nGhost[1];
371: for (PetscInt idx=0; idx<n; ++idx) {
372: const PetscInt eLocalx = pos[idx].i - stag->startGhost[0];
373: const PetscInt eLocaly = pos[idx].j - stag->startGhost[1];
374: const PetscInt eLocalz = pos[idx].k - stag->startGhost[2];
375: const PetscInt eLocal = epl*eLocalz + epr*eLocaly + eLocalx;
377: ix[idx] = eLocal * epe + stag->locationOffsets[pos[idx].loc] + pos[idx].c;
378: }
379: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
380: return 0;
381: }
383: /*@C
384: DMStagVecGetValuesStencil - get vector values using grid indexing
386: Not Collective
388: Input Parameters:
389: + dm - the DMStag object
390: . vec - the vector object
391: . n - the number of values to obtain
392: - pos - locations to obtain values from (as an array of DMStagStencil values)
394: Output Parameter:
395: . val - value at the point
397: Notes:
398: Accepts stencils which refer to global element numbers, but
399: only allows access to entries in the local representation (including ghosts).
401: This approach is not as efficient as getting values directly with DMStagVecGetArray(), which is recommended for matrix free operators.
403: Level: advanced
405: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecSetValuesStencil(), DMStagMatSetValuesStencil(), DMStagVecGetArray()
406: @*/
407: PetscErrorCode DMStagVecGetValuesStencil(DM dm, Vec vec,PetscInt n,const DMStagStencil *pos,PetscScalar *val)
408: {
409: DM_Stag * const stag = (DM_Stag*)dm->data;
410: PetscInt nLocal,idx;
411: PetscInt *ix;
412: PetscScalar const *arr;
416: VecGetLocalSize(vec,&nLocal);
418: PetscMalloc1(n,&ix);
419: DMStagStencilToIndexLocal(dm,dm->dim,n,pos,ix);
420: VecGetArrayRead(vec,&arr);
421: for (idx=0; idx<n; ++idx) val[idx] = arr[ix[idx]];
422: VecRestoreArrayRead(vec,&arr);
423: PetscFree(ix);
424: return 0;
425: }
427: /*@C
428: DMStagVecSetValuesStencil - Set Vec values using global grid indexing
430: Not Collective
432: Input Parameters:
433: + dm - the DMStag object
434: . vec - the Vec
435: . n - the number of values to set
436: . pos - the locations to set values, as an array of DMStagStencil structs
437: . val - the values to set
438: - insertMode - INSERT_VALUES or ADD_VALUES
440: Notes:
441: The vector is expected to be a global vector compatible with the DM (usually obtained by DMGetGlobalVector() or DMCreateGlobalVector()).
443: This approach is not as efficient as setting values directly with DMStagVecGetArray(), which is recommended for matrix-free operators.
444: 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()).
446: Level: advanced
448: .seealso: DMSTAG, DMStagStencil, DMStagStencilLocation, DMStagVecGetValuesStencil(), DMStagMatSetValuesStencil(), DMCreateGlobalVector(), DMGetLocalVector(), DMStagVecGetArray()
449: @*/
450: PetscErrorCode DMStagVecSetValuesStencil(DM dm,Vec vec,PetscInt n,const DMStagStencil *pos,const PetscScalar *val,InsertMode insertMode)
451: {
452: DM_Stag * const stag = (DM_Stag*)dm->data;
453: PetscInt nLocal;
454: PetscInt *ix;
458: VecGetLocalSize(vec,&nLocal);
460: PetscMalloc1(n,&ix);
461: DMStagStencilToIndexLocal(dm,dm->dim,n,pos,ix);
462: VecSetValuesLocal(vec,n,ix,val,insertMode);
463: PetscFree(ix);
464: return 0;
465: }