Actual source code: vsection.c

petsc-3.8.2 2017-11-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: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
 98: {
 99:   PetscScalar    *baseArray;
100:   const PetscInt p = point - s->pStart;

104:   VecGetArray(v, &baseArray);
105:   *values = &baseArray[s->atlasOff[p]];
106:   VecRestoreArray(v, &baseArray);
107:   return(0);
108: }

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

113:   Not collective

115:   Input Parameters:
116: + v - the Vec
117: . s - the organizing PetscSection
118: . point - the point
119: . values - the array of input values
120: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES

122:   Level: developer

124:   Note: This is similar to MatSetValuesStencil(). The Fortran binding is
125: $
126: $   VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
127: $

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

143:   PetscSectionGetConstraintDof(s, point, &cDim);
144:   VecGetArray(v, &baseArray);
145:   array = &baseArray[s->atlasOff[p]];
146:   if (!cDim && doInterior) {
147:     if (orientation >= 0) {
148:       const PetscInt dim = s->atlasDof[p];
149:       PetscInt       i;

151:       if (doInsert) {
152:         for (i = 0; i < dim; ++i) array[i] = values[i];
153:       } else {
154:         for (i = 0; i < dim; ++i) array[i] += values[i];
155:       }
156:     } else {
157:       PetscInt offset = 0;
158:       PetscInt j      = -1, field, i;

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

163:         for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
164:         offset += dim;
165:       }
166:     }
167:   } else if (cDim) {
168:     if (orientation >= 0) {
169:       const PetscInt dim  = s->atlasDof[p];
170:       PetscInt       cInd = 0, i;
171:       const PetscInt *cDof;

173:       PetscSectionGetConstraintIndices(s, point, &cDof);
174:       if (doInsert) {
175:         for (i = 0; i < dim; ++i) {
176:           if ((cInd < cDim) && (i == cDof[cInd])) {
177:             if (doBC) array[i] = values[i]; /* Constrained update */
178:             ++cInd;
179:             continue;
180:           }
181:           if (doInterior) array[i] = values[i]; /* Unconstrained update */
182:         }
183:       } else {
184:         for (i = 0; i < dim; ++i) {
185:           if ((cInd < cDim) && (i == cDof[cInd])) {
186:             if (doBC) array[i] += values[i]; /* Constrained update */
187:             ++cInd;
188:             continue;
189:           }
190:           if (doInterior) array[i] += values[i]; /* Unconstrained update */
191:         }
192:       }
193:     } else {
194:       /* TODO This is broken for add and constrained update */
195:       const PetscInt *cDof;
196:       PetscInt       offset  = 0;
197:       PetscInt       cOffset = 0;
198:       PetscInt       j       = 0, field;

200:       PetscSectionGetConstraintIndices(s, point, &cDof);
201:       for (field = 0; field < s->numFields; ++field) {
202:         const PetscInt dim  = s->field[field]->atlasDof[p];     /* PetscSectionGetFieldDof() */
203:         const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
204:         const PetscInt sDim = dim - tDim;
205:         PetscInt       cInd = 0, i ,k;

207:         for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
208:           if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
209:           if (doInterior) array[j] = values[k];   /* Unconstrained update */
210:         }
211:         offset  += dim;
212:         cOffset += dim - tDim;
213:       }
214:     }
215:   }
216:   VecRestoreArray(v, &baseArray);
217:   return(0);
218: }

220: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
221: {
222:   PetscInt      *subIndices;
223:   PetscInt       Nc, subSize = 0, subOff = 0, p;

227:   PetscSectionGetFieldComponents(section, field, &Nc);
228:   for (p = pStart; p < pEnd; ++p) {
229:     PetscInt gdof, fdof = 0;

231:     PetscSectionGetDof(sectionGlobal, p, &gdof);
232:     if (gdof > 0) {PetscSectionGetFieldDof(section, p, field, &fdof);}
233:     subSize += fdof;
234:   }
235:   PetscMalloc1(subSize, &subIndices);
236:   for (p = pStart; p < pEnd; ++p) {
237:     PetscInt gdof, goff;

239:     PetscSectionGetDof(sectionGlobal, p, &gdof);
240:     if (gdof > 0) {
241:       PetscInt fdof, fc, f2, poff = 0;

243:       PetscSectionGetOffset(sectionGlobal, p, &goff);
244:       /* Can get rid of this loop by storing field information in the global section */
245:       for (f2 = 0; f2 < field; ++f2) {
246:         PetscSectionGetFieldDof(section, p, f2, &fdof);
247:         poff += fdof;
248:       }
249:       PetscSectionGetFieldDof(section, p, field, &fdof);
250:       for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff+poff+fc;
251:     }
252:   }
253:   ISCreateGeneral(PetscObjectComm((PetscObject) v), subSize, subIndices, PETSC_OWN_POINTER, is);
254:   VecGetSubVector(v, *is, subv);
255:   VecSetBlockSize(*subv, Nc);
256:   return(0);
257: }

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

264:   VecRestoreSubVector(v, *is, subv);
265:   ISDestroy(is);
266:   return(0);
267: }

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

272:   Input Parameters:
273: + s    - the local Section
274: . gs   - the global section
275: . x    - the vector
276: - type - one of NORM_1, NORM_2, NORM_INFINITY.

278:   Output Parameter:
279: . val  - the array of norms

281:   Level: intermediate

283: .seealso: VecNorm(), PetscSectionCreate()
284: @*/
285: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
286: {
287:   PetscInt       Nf, f, pStart, pEnd;

295:   PetscSectionGetNumFields(s, &Nf);
296:   if (Nf < 2) {VecNorm(x, type, val);}
297:   else {
298:     PetscSectionGetChart(s, &pStart, &pEnd);
299:     for (f = 0; f < Nf; ++f) {
300:       Vec subv;
301:       IS  is;

303:       PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
304:       VecNorm(subv, type, &val[f]);
305:       PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
306:     }
307:   }
308:   return(0);
309: }