Actual source code: vsection.c

petsc-3.9.0 2018-04-07
Report Typos and Errors
  1: /*
  2:    This file contains routines for section object operations on Vecs
  3: */
  4:  #include <petsc/private/isimpl.h>
  5:  #include <petsc/private/vecimpl.h>

  7: 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: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
 69: {
 70:   PetscBool      isascii;
 71:   PetscInt       f;

 77:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);}
 79:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
 80:   if (isascii) {
 81:     const char *name;

 83:     PetscObjectGetName((PetscObject) v, &name);
 84:     if (s->numFields) {
 85:       PetscViewerASCIIPrintf(viewer, "%s with %D fields\n", name, s->numFields);
 86:       for (f = 0; f < s->numFields; ++f) {
 87:         PetscViewerASCIIPrintf(viewer, "  field %D with %D components\n", f, s->numFieldComponents[f]);
 88:         PetscSectionVecView_ASCII(s->field[f], v, viewer);
 89:       }
 90:     } else {
 91:       PetscViewerASCIIPrintf(viewer, "%s\n", name);
 92:       PetscSectionVecView_ASCII(s, v, viewer);
 93:     }
 94:   }
 95:   return(0);
 96: }

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

101:   Not collective

103:   Input Parameters:
104: + v - the Vec
105: . s - the organizing PetscSection
106: - point - the point

108:   Output Parameter:
109: . values - the array of output values

111:   Level: developer

113: .seealso: PetscSection, PetscSectionCreate(), VecSetValuesSection()
114: @*/
115: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
116: {
117:   PetscScalar    *baseArray;
118:   const PetscInt p = point - s->pStart;

124:   VecGetArray(v, &baseArray);
125:   *values = &baseArray[s->atlasOff[p]];
126:   VecRestoreArray(v, &baseArray);
127:   return(0);
128: }

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

133:   Not collective

135:   Input Parameters:
136: + v - the Vec
137: . s - the organizing PetscSection
138: . point - the point
139: . values - the array of input values
140: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES

142:   Level: developer

144:   Note: This is similar to MatSetValuesStencil(). The Fortran binding is
145: $
146: $   VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
147: $

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

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

173:       if (doInsert) {
174:         for (i = 0; i < dim; ++i) array[i] = values[i];
175:       } else {
176:         for (i = 0; i < dim; ++i) array[i] += values[i];
177:       }
178:     } else {
179:       PetscInt offset = 0;
180:       PetscInt j      = -1, field, i;

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

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

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

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

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

242: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
243: {
244:   PetscInt      *subIndices;
245:   PetscInt       Nc, subSize = 0, subOff = 0, p;

249:   PetscSectionGetFieldComponents(section, field, &Nc);
250:   for (p = pStart; p < pEnd; ++p) {
251:     PetscInt gdof, fdof = 0;

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

261:     PetscSectionGetDof(sectionGlobal, p, &gdof);
262:     if (gdof > 0) {
263:       PetscInt fdof, fc, f2, poff = 0;

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

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

286:   VecRestoreSubVector(v, *is, subv);
287:   ISDestroy(is);
288:   return(0);
289: }

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

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

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

303:   Level: intermediate

305: .seealso: VecNorm(), PetscSectionCreate()
306: @*/
307: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
308: {
309:   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: }