Actual source code: vsection.c
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;
13: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
14: VecGetArray(v, &array);
15: PetscViewerASCIIPushSynchronized(viewer);
16: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
17: for (p = 0; p < s->pEnd - s->pStart; ++p) {
18: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
19: PetscInt b;
21: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT, p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
22: for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
23: PetscScalar v = array[i];
24: #if defined(PETSC_USE_COMPLEX)
25: if (PetscImaginaryPart(v) > 0.0) {
26: PetscViewerASCIISynchronizedPrintf(viewer," %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));
27: } else if (PetscImaginaryPart(v) < 0.0) {
28: PetscViewerASCIISynchronizedPrintf(viewer," %g - %g i", (double)PetscRealPart(v),(double)(-PetscImaginaryPart(v)));
29: } else {
30: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));
31: }
32: #else
33: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);
34: #endif
35: }
36: PetscViewerASCIISynchronizedPrintf(viewer, " constrained");
37: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
38: PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p]+b]);
39: }
40: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
41: } else {
42: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT, p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
43: for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
44: PetscScalar v = array[i];
45: #if defined(PETSC_USE_COMPLEX)
46: if (PetscImaginaryPart(v) > 0.0) {
47: PetscViewerASCIISynchronizedPrintf(viewer," %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));
48: } else if (PetscImaginaryPart(v) < 0.0) {
49: PetscViewerASCIISynchronizedPrintf(viewer," %g - %g i", (double)PetscRealPart(v),(double)(-PetscImaginaryPart(v)));
50: } else {
51: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));
52: }
53: #else
54: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);
55: #endif
56: }
57: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
58: }
59: }
60: PetscViewerFlush(viewer);
61: PetscViewerASCIIPopSynchronized(viewer);
62: VecRestoreArray(v, &array);
63: return 0;
64: }
66: /*@
67: PetscSectionVecView - View a vector, using the section to structure the values
69: Not collective
71: Input Parameters:
72: + s - the organizing PetscSection
73: . v - the Vec
74: - viewer - the PetscViewer
76: Level: developer
78: .seealso: PetscSection, PetscSectionCreate(), VecSetValuesSection()
79: @*/
80: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
81: {
82: PetscBool isascii;
83: PetscInt f;
87: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);
89: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
90: if (isascii) {
91: const char *name;
93: PetscObjectGetName((PetscObject) v, &name);
94: if (s->numFields) {
95: PetscViewerASCIIPrintf(viewer, "%s with %" PetscInt_FMT " fields\n", name, s->numFields);
96: for (f = 0; f < s->numFields; ++f) {
97: PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
98: PetscSectionVecView_ASCII(s->field[f], v, viewer);
99: }
100: } else {
101: PetscViewerASCIIPrintf(viewer, "%s\n", name);
102: PetscSectionVecView_ASCII(s, v, viewer);
103: }
104: }
105: return 0;
106: }
108: /*@C
109: VecGetValuesSection - Gets all the values associated with a given point, according to the section, in the given Vec
111: Not collective
113: Input Parameters:
114: + v - the Vec
115: . s - the organizing PetscSection
116: - point - the point
118: Output Parameter:
119: . values - the array of output values
121: Level: developer
123: .seealso: PetscSection, PetscSectionCreate(), VecSetValuesSection()
124: @*/
125: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
126: {
127: PetscScalar *baseArray;
128: const PetscInt p = point - s->pStart;
132: VecGetArray(v, &baseArray);
133: *values = &baseArray[s->atlasOff[p]];
134: VecRestoreArray(v, &baseArray);
135: return 0;
136: }
138: /*@C
139: VecSetValuesSection - Sets all the values associated with a given point, according to the section, in the given Vec
141: Not collective
143: Input Parameters:
144: + v - the Vec
145: . s - the organizing PetscSection
146: . point - the point
147: . values - the array of input values
148: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES
150: Level: developer
152: Note: This is similar to MatSetValuesStencil(). The Fortran binding is
153: $
154: $ VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
155: $
157: .seealso: PetscSection, PetscSectionCreate(), VecGetValuesSection()
158: @*/
159: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
160: {
161: PetscScalar *baseArray, *array;
162: const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
163: const PetscBool doInterior = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
164: const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
165: const PetscInt p = point - s->pStart;
166: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
167: PetscInt cDim = 0;
171: PetscSectionGetConstraintDof(s, point, &cDim);
172: VecGetArray(v, &baseArray);
173: array = &baseArray[s->atlasOff[p]];
174: if (!cDim && doInterior) {
175: if (orientation >= 0) {
176: const PetscInt dim = s->atlasDof[p];
177: PetscInt i;
179: if (doInsert) {
180: for (i = 0; i < dim; ++i) array[i] = values[i];
181: } else {
182: for (i = 0; i < dim; ++i) array[i] += values[i];
183: }
184: } else {
185: PetscInt offset = 0;
186: PetscInt j = -1, field, i;
188: for (field = 0; field < s->numFields; ++field) {
189: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
191: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
192: offset += dim;
193: }
194: }
195: } else if (cDim) {
196: if (orientation >= 0) {
197: const PetscInt dim = s->atlasDof[p];
198: PetscInt cInd = 0, i;
199: const PetscInt *cDof;
201: PetscSectionGetConstraintIndices(s, point, &cDof);
202: if (doInsert) {
203: for (i = 0; i < dim; ++i) {
204: if ((cInd < cDim) && (i == cDof[cInd])) {
205: if (doBC) array[i] = values[i]; /* Constrained update */
206: ++cInd;
207: continue;
208: }
209: if (doInterior) array[i] = values[i]; /* Unconstrained update */
210: }
211: } else {
212: for (i = 0; i < dim; ++i) {
213: if ((cInd < cDim) && (i == cDof[cInd])) {
214: if (doBC) array[i] += values[i]; /* Constrained update */
215: ++cInd;
216: continue;
217: }
218: if (doInterior) array[i] += values[i]; /* Unconstrained update */
219: }
220: }
221: } else {
222: /* TODO This is broken for add and constrained update */
223: const PetscInt *cDof;
224: PetscInt offset = 0;
225: PetscInt cOffset = 0;
226: PetscInt j = 0, field;
228: PetscSectionGetConstraintIndices(s, point, &cDof);
229: for (field = 0; field < s->numFields; ++field) {
230: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
231: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
232: const PetscInt sDim = dim - tDim;
233: PetscInt cInd = 0, i ,k;
235: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
236: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
237: if (doInterior) array[j] = values[k]; /* Unconstrained update */
238: }
239: offset += dim;
240: cOffset += dim - tDim;
241: }
242: }
243: }
244: VecRestoreArray(v, &baseArray);
245: return 0;
246: }
248: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
249: {
250: PetscInt *subIndices;
251: PetscInt Nc, subSize = 0, subOff = 0, p;
253: PetscSectionGetFieldComponents(section, field, &Nc);
254: for (p = pStart; p < pEnd; ++p) {
255: PetscInt gdof, fdof = 0;
257: PetscSectionGetDof(sectionGlobal, p, &gdof);
258: if (gdof > 0) PetscSectionGetFieldDof(section, p, field, &fdof);
259: subSize += fdof;
260: }
261: PetscMalloc1(subSize, &subIndices);
262: for (p = pStart; p < pEnd; ++p) {
263: PetscInt gdof, goff;
265: PetscSectionGetDof(sectionGlobal, p, &gdof);
266: if (gdof > 0) {
267: PetscInt fdof, fc, f2, poff = 0;
269: PetscSectionGetOffset(sectionGlobal, p, &goff);
270: /* Can get rid of this loop by storing field information in the global section */
271: for (f2 = 0; f2 < field; ++f2) {
272: PetscSectionGetFieldDof(section, p, f2, &fdof);
273: poff += fdof;
274: }
275: PetscSectionGetFieldDof(section, p, field, &fdof);
276: for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff+poff+fc;
277: }
278: }
279: ISCreateGeneral(PetscObjectComm((PetscObject) v), subSize, subIndices, PETSC_OWN_POINTER, is);
280: VecGetSubVector(v, *is, subv);
281: VecSetBlockSize(*subv, Nc);
282: return 0;
283: }
285: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
286: {
287: VecRestoreSubVector(v, *is, subv);
288: ISDestroy(is);
289: return 0;
290: }
292: /*@C
293: PetscSectionVecNorm - Computes the vector norm, separated into field components.
295: Input Parameters:
296: + s - the local Section
297: . gs - the global section
298: . x - the vector
299: - type - one of NORM_1, NORM_2, NORM_INFINITY.
301: Output Parameter:
302: . val - the array of norms
304: Level: intermediate
306: .seealso: VecNorm(), PetscSectionCreate()
307: @*/
308: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
309: {
310: PetscInt Nf, f, pStart, pEnd;
316: PetscSectionGetNumFields(s, &Nf);
317: if (Nf < 2) VecNorm(x, type, val);
318: else {
319: PetscSectionGetChart(s, &pStart, &pEnd);
320: for (f = 0; f < Nf; ++f) {
321: Vec subv;
322: IS is;
324: PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
325: VecNorm(subv, type, &val[f]);
326: PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
327: }
328: }
329: return 0;
330: }