Actual source code: vsection.c

petsc-3.7.3 2016-07-24
Report Typos and Errors
  1: /*
  2:    This file contains routines for section object operations on Vecs
  3: */
  4: #include <petsc/private/isimpl.h>   /*I  "petscvec.h"   I*/
  5: #include <petsc/private/vecimpl.h>   /*I  "petscvec.h"   I*/

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

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

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

 72: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
 73: {
 74:   PetscBool      isascii;
 75:   PetscInt       f;

 79:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);}
 82:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
 83:   if (isascii) {
 84:     const char *name;

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

103: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
104: {
105:   PetscScalar    *baseArray;
106:   const PetscInt p = point - s->pStart;

110:   VecGetArray(v, &baseArray);
111:   *values = &baseArray[s->atlasOff[p]];
112:   VecRestoreArray(v, &baseArray);
113:   return(0);
114: }

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

121:   Not collective

123:   Input Parameters:
124: + v - the Vec
125: . s - the organizing PetscSection
126: . point - the point
127: . values - the array of input values
128: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES

130:   Level: developer

132:   Note: This is similar to MatSetValuesStencil(). The Fortran binding is
133: $
134: $   VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
135: $

137: .seealso: PetscSection, PetscSectionCreate()
138: @*/
139: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
140: {
141:   PetscScalar     *baseArray, *array;
142:   const PetscBool doInsert    = mode == INSERT_VALUES     || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES                          ? PETSC_TRUE : PETSC_FALSE;
143:   const PetscBool doInterior  = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES    || mode == INSERT_VALUES    || mode == ADD_VALUES    ? PETSC_TRUE : PETSC_FALSE;
144:   const PetscBool doBC        = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES    || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
145:   const PetscInt  p           = point - s->pStart;
146:   const PetscInt  orientation = 0; /* Needs to be included for use in closure operations */
147:   PetscInt        cDim        = 0;
148:   PetscErrorCode  ierr;

151:   PetscSectionGetConstraintDof(s, point, &cDim);
152:   VecGetArray(v, &baseArray);
153:   array = &baseArray[s->atlasOff[p]];
154:   if (!cDim && doInterior) {
155:     if (orientation >= 0) {
156:       const PetscInt dim = s->atlasDof[p];
157:       PetscInt       i;

159:       if (doInsert) {
160:         for (i = 0; i < dim; ++i) array[i] = values[i];
161:       } else {
162:         for (i = 0; i < dim; ++i) array[i] += values[i];
163:       }
164:     } else {
165:       PetscInt offset = 0;
166:       PetscInt j      = -1, field, i;

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

171:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
172:         offset += dim;
173:       }
174:     }
175:   } else if (cDim) {
176:     if (orientation >= 0) {
177:       const PetscInt dim  = s->atlasDof[p];
178:       PetscInt       cInd = 0, i;
179:       const PetscInt *cDof;

181:       PetscSectionGetConstraintIndices(s, point, &cDof);
182:       if (doInsert) {
183:         for (i = 0; i < dim; ++i) {
184:           if ((cInd < cDim) && (i == cDof[cInd])) {
185:             if (doBC) array[i] = values[i]; /* Constrained update */
186:             ++cInd;
187:             continue;
188:           }
189:           if (doInterior) array[i] = values[i]; /* Unconstrained update */
190:         }
191:       } else {
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:       }
201:     } else {
202:       /* TODO This is broken for add and constrained update */
203:       const PetscInt *cDof;
204:       PetscInt       offset  = 0;
205:       PetscInt       cOffset = 0;
206:       PetscInt       j       = 0, field;

208:       PetscSectionGetConstraintIndices(s, point, &cDof);
209:       for (field = 0; field < s->numFields; ++field) {
210:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
211:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
212:         const PetscInt sDim = dim - tDim;
213:         PetscInt       cInd = 0, i ,k;

215:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
216:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
217:           if (doInterior) array[j] = values[k];   /* Unconstrained update */
218:         }
219:         offset  += dim;
220:         cOffset += dim - tDim;
221:       }
222:     }
223:   }
224:   VecRestoreArray(v, &baseArray);
225:   return(0);
226: }

230: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
231: {
232:   PetscInt      *subIndices;
233:   PetscInt       Nc, subSize = 0, subOff = 0, p;

237:   PetscSectionGetFieldComponents(section, field, &Nc);
238:   for (p = pStart; p < pEnd; ++p) {
239:     PetscInt gdof, fdof = 0;

241:     PetscSectionGetDof(sectionGlobal, p, &gdof);
242:     if (gdof > 0) {PetscSectionGetFieldDof(section, p, field, &fdof);}
243:     subSize += fdof;
244:   }
245:   PetscMalloc1(subSize, &subIndices);
246:   for (p = pStart; p < pEnd; ++p) {
247:     PetscInt gdof, goff;

249:     PetscSectionGetDof(sectionGlobal, p, &gdof);
250:     if (gdof > 0) {
251:       PetscInt fdof, fc, f2, poff = 0;

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

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

276:   VecRestoreSubVector(v, *is, subv);
277:   ISDestroy(is);
278:   return(0);
279: }

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

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

292:   Output Parameter:
293: . val  - the array of norms

295:   Level: intermediate

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

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

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