Actual source code: vsection.c

petsc-3.8.3 2017-12-09
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;

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

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

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

100:   Not collective

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

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

110:   Level: developer

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

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

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

130:   Not collective

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

139:   Level: developer

141:   Note: This is similar to MatSetValuesStencil(). The Fortran binding is
142: $
143: $   VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
144: $

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

160:   PetscSectionGetConstraintDof(s, point, &cDim);
161:   VecGetArray(v, &baseArray);
162:   array = &baseArray[s->atlasOff[p]];
163:   if (!cDim && doInterior) {
164:     if (orientation >= 0) {
165:       const PetscInt dim = s->atlasDof[p];
166:       PetscInt       i;

168:       if (doInsert) {
169:         for (i = 0; i < dim; ++i) array[i] = values[i];
170:       } else {
171:         for (i = 0; i < dim; ++i) array[i] += values[i];
172:       }
173:     } else {
174:       PetscInt offset = 0;
175:       PetscInt j      = -1, field, i;

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

180:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
181:         offset += dim;
182:       }
183:     }
184:   } else if (cDim) {
185:     if (orientation >= 0) {
186:       const PetscInt dim  = s->atlasDof[p];
187:       PetscInt       cInd = 0, i;
188:       const PetscInt *cDof;

190:       PetscSectionGetConstraintIndices(s, point, &cDof);
191:       if (doInsert) {
192:         for (i = 0; i < dim; ++i) {
193:           if ((cInd < cDim) && (i == cDof[cInd])) {
194:             if (doBC) array[i] = values[i]; /* Constrained update */
195:             ++cInd;
196:             continue;
197:           }
198:           if (doInterior) array[i] = values[i]; /* Unconstrained update */
199:         }
200:       } else {
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:       }
210:     } else {
211:       /* TODO This is broken for add and constrained update */
212:       const PetscInt *cDof;
213:       PetscInt       offset  = 0;
214:       PetscInt       cOffset = 0;
215:       PetscInt       j       = 0, field;

217:       PetscSectionGetConstraintIndices(s, point, &cDof);
218:       for (field = 0; field < s->numFields; ++field) {
219:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
220:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
221:         const PetscInt sDim = dim - tDim;
222:         PetscInt       cInd = 0, i ,k;

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

237: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
238: {
239:   PetscInt      *subIndices;
240:   PetscInt       Nc, subSize = 0, subOff = 0, p;

244:   PetscSectionGetFieldComponents(section, field, &Nc);
245:   for (p = pStart; p < pEnd; ++p) {
246:     PetscInt gdof, fdof = 0;

248:     PetscSectionGetDof(sectionGlobal, p, &gdof);
249:     if (gdof > 0) {PetscSectionGetFieldDof(section, p, field, &fdof);}
250:     subSize += fdof;
251:   }
252:   PetscMalloc1(subSize, &subIndices);
253:   for (p = pStart; p < pEnd; ++p) {
254:     PetscInt gdof, goff;

256:     PetscSectionGetDof(sectionGlobal, p, &gdof);
257:     if (gdof > 0) {
258:       PetscInt fdof, fc, f2, poff = 0;

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

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

281:   VecRestoreSubVector(v, *is, subv);
282:   ISDestroy(is);
283:   return(0);
284: }

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

289:   Input Parameters:
290: + s    - the local Section
291: . gs   - the global section
292: . x    - the vector
293: - type - one of NORM_1, NORM_2, NORM_INFINITY.

295:   Output Parameter:
296: . val  - the array of norms

298:   Level: intermediate

300: .seealso: VecNorm(), PetscSectionCreate()
301: @*/
302: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
303: {
304:   PetscInt       Nf, f, pStart, pEnd;

312:   PetscSectionGetNumFields(s, &Nf);
313:   if (Nf < 2) {VecNorm(x, type, val);}
314:   else {
315:     PetscSectionGetChart(s, &pStart, &pEnd);
316:     for (f = 0; f < Nf; ++f) {
317:       Vec subv;
318:       IS  is;

320:       PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
321:       VecNorm(subv, type, &val[f]);
322:       PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
323:     }
324:   }
325:   return(0);
326: }