Actual source code: vsection.c

  1: /*
  2:    This file contains routines for section object operations on Vecs
  3: */
  4: #include <petsc/private/sectionimpl.h>
  5: #include <petsc/private/vecimpl.h>

  7: static PetscErrorCode PetscSectionVecView_ASCII(PetscSection s, Vec v, PetscViewer viewer)
  8: {
  9:   PetscScalar *array;
 10:   PetscInt     p, i;
 11:   PetscMPIInt  rank;

 13:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
 14:   VecGetArray(v, &array);
 15:   PetscViewerASCIIPushSynchronized(viewer);
 16:   PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
 17:   for (p = 0; p < s->pEnd - s->pStart; ++p) {
 18:     if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
 19:       PetscInt b;

 21:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
 22:       for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) {
 23:         PetscScalar v = array[i];
 24: #if defined(PETSC_USE_COMPLEX)
 25:         if (PetscImaginaryPart(v) > 0.0) {
 26:           PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));
 27:         } else if (PetscImaginaryPart(v) < 0.0) {
 28:           PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPart(v), (double)(-PetscImaginaryPart(v)));
 29:         } else {
 30:           PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));
 31:         }
 32: #else
 33:         PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);
 34: #endif
 35:       }
 36:       PetscViewerASCIISynchronizedPrintf(viewer, " constrained");
 37:       for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]);
 38:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
 39:     } else {
 40:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
 41:       for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) {
 42:         PetscScalar v = array[i];
 43: #if defined(PETSC_USE_COMPLEX)
 44:         if (PetscImaginaryPart(v) > 0.0) {
 45:           PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));
 46:         } else if (PetscImaginaryPart(v) < 0.0) {
 47:           PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPart(v), (double)(-PetscImaginaryPart(v)));
 48:         } else {
 49:           PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));
 50:         }
 51: #else
 52:         PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);
 53: #endif
 54:       }
 55:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
 56:     }
 57:   }
 58:   PetscViewerFlush(viewer);
 59:   PetscViewerASCIIPopSynchronized(viewer);
 60:   VecRestoreArray(v, &array);
 61:   return 0;
 62: }

 64: /*@
 65:   PetscSectionVecView - View a vector, using the section to structure the values

 67:   Not collective

 69:   Input Parameters:
 70: + s      - the organizing PetscSection
 71: . v      - the Vec
 72: - viewer - the PetscViewer

 74:   Level: developer

 76: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecSetValuesSection()`
 77: @*/
 78: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
 79: {
 80:   PetscBool isascii;
 81:   PetscInt  f;

 85:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);
 87:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
 88:   if (isascii) {
 89:     const char *name;

 91:     PetscObjectGetName((PetscObject)v, &name);
 92:     if (s->numFields) {
 93:       PetscViewerASCIIPrintf(viewer, "%s with %" PetscInt_FMT " fields\n", name, s->numFields);
 94:       for (f = 0; f < s->numFields; ++f) {
 95:         PetscViewerASCIIPrintf(viewer, "  field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
 96:         PetscSectionVecView_ASCII(s->field[f], v, viewer);
 97:       }
 98:     } else {
 99:       PetscViewerASCIIPrintf(viewer, "%s\n", name);
100:       PetscSectionVecView_ASCII(s, v, viewer);
101:     }
102:   }
103:   return 0;
104: }

106: /*@C
107:   VecGetValuesSection - Gets all the values associated with a given point, according to the section, in the given Vec

109:   Not collective

111:   Input Parameters:
112: + v - the Vec
113: . s - the organizing PetscSection
114: - point - the point

116:   Output Parameter:
117: . values - the array of output values

119:   Level: developer

121: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecSetValuesSection()`
122: @*/
123: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
124: {
125:   PetscScalar   *baseArray;
126:   const PetscInt p = point - s->pStart;

130:   VecGetArray(v, &baseArray);
131:   *values = &baseArray[s->atlasOff[p]];
132:   VecRestoreArray(v, &baseArray);
133:   return 0;
134: }

136: /*@C
137:   VecSetValuesSection - Sets all the values associated with a given point, according to the section, in the given Vec

139:   Not collective

141:   Input Parameters:
142: + v - the Vec
143: . s - the organizing PetscSection
144: . point - the point
145: . values - the array of input values
146: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES

148:   Level: developer

150:   Note: This is similar to MatSetValuesStencil(). The Fortran binding is
151: $
152: $   VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
153: $

155: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecGetValuesSection()`
156: @*/
157: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
158: {
159:   PetscScalar    *baseArray, *array;
160:   const PetscBool doInsert    = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
161:   const PetscBool doInterior  = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
162:   const PetscBool doBC        = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
163:   const PetscInt  p           = point - s->pStart;
164:   const PetscInt  orientation = 0; /* Needs to be included for use in closure operations */
165:   PetscInt        cDim        = 0;

169:   PetscSectionGetConstraintDof(s, point, &cDim);
170:   VecGetArray(v, &baseArray);
171:   array = &baseArray[s->atlasOff[p]];
172:   if (!cDim && doInterior) {
173:     if (orientation >= 0) {
174:       const PetscInt dim = s->atlasDof[p];
175:       PetscInt       i;

177:       if (doInsert) {
178:         for (i = 0; i < dim; ++i) array[i] = values[i];
179:       } else {
180:         for (i = 0; i < dim; ++i) array[i] += values[i];
181:       }
182:     } else {
183:       PetscInt offset = 0;
184:       PetscInt j      = -1, field, i;

186:       for (field = 0; field < s->numFields; ++field) {
187:         const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */

189:         for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
190:         offset += dim;
191:       }
192:     }
193:   } else if (cDim) {
194:     if (orientation >= 0) {
195:       const PetscInt  dim  = s->atlasDof[p];
196:       PetscInt        cInd = 0, i;
197:       const PetscInt *cDof;

199:       PetscSectionGetConstraintIndices(s, point, &cDof);
200:       if (doInsert) {
201:         for (i = 0; i < dim; ++i) {
202:           if ((cInd < cDim) && (i == cDof[cInd])) {
203:             if (doBC) array[i] = values[i]; /* Constrained update */
204:             ++cInd;
205:             continue;
206:           }
207:           if (doInterior) array[i] = values[i]; /* Unconstrained update */
208:         }
209:       } else {
210:         for (i = 0; i < dim; ++i) {
211:           if ((cInd < cDim) && (i == cDof[cInd])) {
212:             if (doBC) array[i] += values[i]; /* Constrained update */
213:             ++cInd;
214:             continue;
215:           }
216:           if (doInterior) array[i] += values[i]; /* Unconstrained update */
217:         }
218:       }
219:     } else {
220:       /* TODO This is broken for add and constrained update */
221:       const PetscInt *cDof;
222:       PetscInt        offset  = 0;
223:       PetscInt        cOffset = 0;
224:       PetscInt        j       = 0, field;

226:       PetscSectionGetConstraintIndices(s, point, &cDof);
227:       for (field = 0; field < s->numFields; ++field) {
228:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
229:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
230:         const PetscInt sDim = dim - tDim;
231:         PetscInt       cInd = 0, i, k;

233:         for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
234:           if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
235:             ++cInd;
236:             continue;
237:           }
238:           if (doInterior) array[j] = values[k]; /* Unconstrained update */
239:         }
240:         offset += dim;
241:         cOffset += dim - tDim;
242:       }
243:     }
244:   }
245:   VecRestoreArray(v, &baseArray);
246:   return 0;
247: }

249: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
250: {
251:   PetscInt *subIndices;
252:   PetscInt  Nc, subSize = 0, subOff = 0, p;

254:   PetscSectionGetFieldComponents(section, field, &Nc);
255:   for (p = pStart; p < pEnd; ++p) {
256:     PetscInt gdof, fdof = 0;

258:     PetscSectionGetDof(sectionGlobal, p, &gdof);
259:     if (gdof > 0) PetscSectionGetFieldDof(section, p, field, &fdof);
260:     subSize += fdof;
261:   }
262:   PetscMalloc1(subSize, &subIndices);
263:   for (p = pStart; p < pEnd; ++p) {
264:     PetscInt gdof, goff;

266:     PetscSectionGetDof(sectionGlobal, p, &gdof);
267:     if (gdof > 0) {
268:       PetscInt fdof, fc, f2, poff = 0;

270:       PetscSectionGetOffset(sectionGlobal, p, &goff);
271:       /* Can get rid of this loop by storing field information in the global section */
272:       for (f2 = 0; f2 < field; ++f2) {
273:         PetscSectionGetFieldDof(section, p, f2, &fdof);
274:         poff += fdof;
275:       }
276:       PetscSectionGetFieldDof(section, p, field, &fdof);
277:       for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
278:     }
279:   }
280:   ISCreateGeneral(PetscObjectComm((PetscObject)v), subSize, subIndices, PETSC_OWN_POINTER, is);
281:   VecGetSubVector(v, *is, subv);
282:   VecSetBlockSize(*subv, Nc);
283:   return 0;
284: }

286: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
287: {
288:   VecRestoreSubVector(v, *is, subv);
289:   ISDestroy(is);
290:   return 0;
291: }

293: /*@C
294:   PetscSectionVecNorm - Computes the vector norm, separated into field components.

296:   Input Parameters:
297: + s    - the local Section
298: . gs   - the global section
299: . x    - the vector
300: - type - one of NORM_1, NORM_2, NORM_INFINITY.

302:   Output Parameter:
303: . val  - the array of norms

305:   Level: intermediate

307: .seealso: `VecNorm()`, `PetscSectionCreate()`
308: @*/
309: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
310: {
311:   PetscInt Nf, f, pStart, pEnd;

317:   PetscSectionGetNumFields(s, &Nf);
318:   if (Nf < 2) VecNorm(x, type, val);
319:   else {
320:     PetscSectionGetChart(s, &pStart, &pEnd);
321:     for (f = 0; f < Nf; ++f) {
322:       Vec subv;
323:       IS  is;

325:       PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
326:       VecNorm(subv, type, &val[f]);
327:       PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
328:     }
329:   }
330:   return 0;
331: }