Actual source code: vsection.c

petsc-master 2019-09-18
Report Typos and Errors
  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;

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

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

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

 71:   Not collective

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

 78:   Level: developer

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

 91:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);}
 93:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
 94:   if (isascii) {
 95:     const char *name;

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

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

115:   Not collective

117:   Input Parameters:
118: + v - the Vec
119: . s - the organizing PetscSection
120: - point - the point

122:   Output Parameter:
123: . values - the array of output values

125:   Level: developer

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

138:   VecGetArray(v, &baseArray);
139:   *values = &baseArray[s->atlasOff[p]];
140:   VecRestoreArray(v, &baseArray);
141:   return(0);
142: }

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

147:   Not collective

149:   Input Parameters:
150: + v - the Vec
151: . s - the organizing PetscSection
152: . point - the point
153: . values - the array of input values
154: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES

156:   Level: developer

158:   Note: This is similar to MatSetValuesStencil(). The Fortran binding is
159: $
160: $   VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
161: $

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

179:   PetscSectionGetConstraintDof(s, point, &cDim);
180:   VecGetArray(v, &baseArray);
181:   array = &baseArray[s->atlasOff[p]];
182:   if (!cDim && doInterior) {
183:     if (orientation >= 0) {
184:       const PetscInt dim = s->atlasDof[p];
185:       PetscInt       i;

187:       if (doInsert) {
188:         for (i = 0; i < dim; ++i) array[i] = values[i];
189:       } else {
190:         for (i = 0; i < dim; ++i) array[i] += values[i];
191:       }
192:     } else {
193:       PetscInt offset = 0;
194:       PetscInt j      = -1, field, i;

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

199:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
200:         offset += dim;
201:       }
202:     }
203:   } else if (cDim) {
204:     if (orientation >= 0) {
205:       const PetscInt dim  = s->atlasDof[p];
206:       PetscInt       cInd = 0, i;
207:       const PetscInt *cDof;

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

236:       PetscSectionGetConstraintIndices(s, point, &cDof);
237:       for (field = 0; field < s->numFields; ++field) {
238:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
239:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
240:         const PetscInt sDim = dim - tDim;
241:         PetscInt       cInd = 0, i ,k;

243:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
244:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
245:           if (doInterior) array[j] = values[k];   /* Unconstrained update */
246:         }
247:         offset  += dim;
248:         cOffset += dim - tDim;
249:       }
250:     }
251:   }
252:   VecRestoreArray(v, &baseArray);
253:   return(0);
254: }

256: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
257: {
258:   PetscInt      *subIndices;
259:   PetscInt       Nc, subSize = 0, subOff = 0, p;

263:   PetscSectionGetFieldComponents(section, field, &Nc);
264:   for (p = pStart; p < pEnd; ++p) {
265:     PetscInt gdof, fdof = 0;

267:     PetscSectionGetDof(sectionGlobal, p, &gdof);
268:     if (gdof > 0) {PetscSectionGetFieldDof(section, p, field, &fdof);}
269:     subSize += fdof;
270:   }
271:   PetscMalloc1(subSize, &subIndices);
272:   for (p = pStart; p < pEnd; ++p) {
273:     PetscInt gdof, goff;

275:     PetscSectionGetDof(sectionGlobal, p, &gdof);
276:     if (gdof > 0) {
277:       PetscInt fdof, fc, f2, poff = 0;

279:       PetscSectionGetOffset(sectionGlobal, p, &goff);
280:       /* Can get rid of this loop by storing field information in the global section */
281:       for (f2 = 0; f2 < field; ++f2) {
282:         PetscSectionGetFieldDof(section, p, f2, &fdof);
283:         poff += fdof;
284:       }
285:       PetscSectionGetFieldDof(section, p, field, &fdof);
286:       for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff+poff+fc;
287:     }
288:   }
289:   ISCreateGeneral(PetscObjectComm((PetscObject) v), subSize, subIndices, PETSC_OWN_POINTER, is);
290:   VecGetSubVector(v, *is, subv);
291:   VecSetBlockSize(*subv, Nc);
292:   return(0);
293: }

295: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
296: {

300:   VecRestoreSubVector(v, *is, subv);
301:   ISDestroy(is);
302:   return(0);
303: }

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

308:   Input Parameters:
309: + s    - the local Section
310: . gs   - the global section
311: . x    - the vector
312: - type - one of NORM_1, NORM_2, NORM_INFINITY.

314:   Output Parameter:
315: . val  - the array of norms

317:   Level: intermediate

319: .seealso: VecNorm(), PetscSectionCreate()
320: @*/
321: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
322: {
323:   PetscInt       Nf, f, pStart, pEnd;

331:   PetscSectionGetNumFields(s, &Nf);
332:   if (Nf < 2) {VecNorm(x, type, val);}
333:   else {
334:     PetscSectionGetChart(s, &pStart, &pEnd);
335:     for (f = 0; f < Nf; ++f) {
336:       Vec subv;
337:       IS  is;

339:       PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
340:       VecNorm(subv, type, &val[f]);
341:       PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
342:     }
343:   }
344:   return(0);
345: }