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) {
 38:         PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p]+b]);
 39:       }
 40:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
 41:     } else {
 42:       PetscViewerASCIISynchronizedPrintf(viewer, "  (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT, p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
 43:       for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
 44:         PetscScalar v = array[i];
 45: #if defined(PETSC_USE_COMPLEX)
 46:         if (PetscImaginaryPart(v) > 0.0) {
 47:           PetscViewerASCIISynchronizedPrintf(viewer," %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));
 48:         } else if (PetscImaginaryPart(v) < 0.0) {
 49:           PetscViewerASCIISynchronizedPrintf(viewer," %g - %g i", (double)PetscRealPart(v),(double)(-PetscImaginaryPart(v)));
 50:         } else {
 51:           PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));
 52:         }
 53: #else
 54:         PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);
 55: #endif
 56:       }
 57:       PetscViewerASCIISynchronizedPrintf(viewer, "\n");
 58:     }
 59:   }
 60:   PetscViewerFlush(viewer);
 61:   PetscViewerASCIIPopSynchronized(viewer);
 62:   VecRestoreArray(v, &array);
 63:   return 0;
 64: }

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

 69:   Not collective

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

 76:   Level: developer

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

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

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

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

111:   Not collective

113:   Input Parameters:
114: + v - the Vec
115: . s - the organizing PetscSection
116: - point - the point

118:   Output Parameter:
119: . values - the array of output values

121:   Level: developer

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

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

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

141:   Not collective

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

150:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

304:   Level: intermediate

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

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

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