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) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]);
38: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
39: } else {
40: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
41: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) {
42: PetscScalar v = array[i];
43: #if defined(PETSC_USE_COMPLEX)
44: if (PetscImaginaryPart(v) > 0.0) {
45: PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));
46: } else if (PetscImaginaryPart(v) < 0.0) {
47: PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPart(v), (double)(-PetscImaginaryPart(v)));
48: } else {
49: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));
50: }
51: #else
52: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);
53: #endif
54: }
55: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
56: }
57: }
58: PetscViewerFlush(viewer);
59: PetscViewerASCIIPopSynchronized(viewer);
60: VecRestoreArray(v, &array);
61: return 0;
62: }
64: /*@
65: PetscSectionVecView - View a vector, using the section to structure the values
67: Not collective
69: Input Parameters:
70: + s - the organizing PetscSection
71: . v - the Vec
72: - viewer - the PetscViewer
74: Level: developer
76: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecSetValuesSection()`
77: @*/
78: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
79: {
80: PetscBool isascii;
81: PetscInt f;
85: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);
87: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
88: if (isascii) {
89: const char *name;
91: PetscObjectGetName((PetscObject)v, &name);
92: if (s->numFields) {
93: PetscViewerASCIIPrintf(viewer, "%s with %" PetscInt_FMT " fields\n", name, s->numFields);
94: for (f = 0; f < s->numFields; ++f) {
95: PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
96: PetscSectionVecView_ASCII(s->field[f], v, viewer);
97: }
98: } else {
99: PetscViewerASCIIPrintf(viewer, "%s\n", name);
100: PetscSectionVecView_ASCII(s, v, viewer);
101: }
102: }
103: return 0;
104: }
106: /*@C
107: VecGetValuesSection - Gets all the values associated with a given point, according to the section, in the given Vec
109: Not collective
111: Input Parameters:
112: + v - the Vec
113: . s - the organizing PetscSection
114: - point - the point
116: Output Parameter:
117: . values - the array of output values
119: Level: developer
121: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecSetValuesSection()`
122: @*/
123: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
124: {
125: PetscScalar *baseArray;
126: const PetscInt p = point - s->pStart;
130: VecGetArray(v, &baseArray);
131: *values = &baseArray[s->atlasOff[p]];
132: VecRestoreArray(v, &baseArray);
133: return 0;
134: }
136: /*@C
137: VecSetValuesSection - Sets all the values associated with a given point, according to the section, in the given Vec
139: Not collective
141: Input Parameters:
142: + v - the Vec
143: . s - the organizing PetscSection
144: . point - the point
145: . values - the array of input values
146: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES
148: Level: developer
150: Note: This is similar to MatSetValuesStencil(). The Fortran binding is
151: $
152: $ VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
153: $
155: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecGetValuesSection()`
156: @*/
157: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
158: {
159: PetscScalar *baseArray, *array;
160: const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
161: const PetscBool doInterior = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
162: const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
163: const PetscInt p = point - s->pStart;
164: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
165: PetscInt cDim = 0;
169: PetscSectionGetConstraintDof(s, point, &cDim);
170: VecGetArray(v, &baseArray);
171: array = &baseArray[s->atlasOff[p]];
172: if (!cDim && doInterior) {
173: if (orientation >= 0) {
174: const PetscInt dim = s->atlasDof[p];
175: PetscInt i;
177: if (doInsert) {
178: for (i = 0; i < dim; ++i) array[i] = values[i];
179: } else {
180: for (i = 0; i < dim; ++i) array[i] += values[i];
181: }
182: } else {
183: PetscInt offset = 0;
184: PetscInt j = -1, field, i;
186: for (field = 0; field < s->numFields; ++field) {
187: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
189: for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
190: offset += dim;
191: }
192: }
193: } else if (cDim) {
194: if (orientation >= 0) {
195: const PetscInt dim = s->atlasDof[p];
196: PetscInt cInd = 0, i;
197: const PetscInt *cDof;
199: PetscSectionGetConstraintIndices(s, point, &cDof);
200: if (doInsert) {
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: } else {
210: for (i = 0; i < dim; ++i) {
211: if ((cInd < cDim) && (i == cDof[cInd])) {
212: if (doBC) array[i] += values[i]; /* Constrained update */
213: ++cInd;
214: continue;
215: }
216: if (doInterior) array[i] += values[i]; /* Unconstrained update */
217: }
218: }
219: } else {
220: /* TODO This is broken for add and constrained update */
221: const PetscInt *cDof;
222: PetscInt offset = 0;
223: PetscInt cOffset = 0;
224: PetscInt j = 0, field;
226: PetscSectionGetConstraintIndices(s, point, &cDof);
227: for (field = 0; field < s->numFields; ++field) {
228: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
229: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
230: const PetscInt sDim = dim - tDim;
231: PetscInt cInd = 0, i, k;
233: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
234: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
235: ++cInd;
236: continue;
237: }
238: if (doInterior) array[j] = values[k]; /* Unconstrained update */
239: }
240: offset += dim;
241: cOffset += dim - tDim;
242: }
243: }
244: }
245: VecRestoreArray(v, &baseArray);
246: return 0;
247: }
249: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
250: {
251: PetscInt *subIndices;
252: PetscInt Nc, subSize = 0, subOff = 0, p;
254: PetscSectionGetFieldComponents(section, field, &Nc);
255: for (p = pStart; p < pEnd; ++p) {
256: PetscInt gdof, fdof = 0;
258: PetscSectionGetDof(sectionGlobal, p, &gdof);
259: if (gdof > 0) PetscSectionGetFieldDof(section, p, field, &fdof);
260: subSize += fdof;
261: }
262: PetscMalloc1(subSize, &subIndices);
263: for (p = pStart; p < pEnd; ++p) {
264: PetscInt gdof, goff;
266: PetscSectionGetDof(sectionGlobal, p, &gdof);
267: if (gdof > 0) {
268: PetscInt fdof, fc, f2, poff = 0;
270: PetscSectionGetOffset(sectionGlobal, p, &goff);
271: /* Can get rid of this loop by storing field information in the global section */
272: for (f2 = 0; f2 < field; ++f2) {
273: PetscSectionGetFieldDof(section, p, f2, &fdof);
274: poff += fdof;
275: }
276: PetscSectionGetFieldDof(section, p, field, &fdof);
277: for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
278: }
279: }
280: ISCreateGeneral(PetscObjectComm((PetscObject)v), subSize, subIndices, PETSC_OWN_POINTER, is);
281: VecGetSubVector(v, *is, subv);
282: VecSetBlockSize(*subv, Nc);
283: return 0;
284: }
286: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
287: {
288: VecRestoreSubVector(v, *is, subv);
289: ISDestroy(is);
290: return 0;
291: }
293: /*@C
294: PetscSectionVecNorm - Computes the vector norm, separated into field components.
296: Input Parameters:
297: + s - the local Section
298: . gs - the global section
299: . x - the vector
300: - type - one of NORM_1, NORM_2, NORM_INFINITY.
302: Output Parameter:
303: . val - the array of norms
305: Level: intermediate
307: .seealso: `VecNorm()`, `PetscSectionCreate()`
308: @*/
309: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
310: {
311: 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: }