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: PetscFunctionBegin;
14: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
15: PetscCall(VecGetArray(v, &array));
16: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
17: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
18: for (p = 0; p < s->pEnd - s->pStart; ++p) {
19: if (s->bc && (s->bc->atlasDof[p] > 0)) {
20: PetscInt b;
22: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
23: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) {
24: PetscScalar v = array[i];
25: #if defined(PETSC_USE_COMPLEX)
26: if (PetscImaginaryPart(v) > 0.0) {
27: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v)));
28: } else if (PetscImaginaryPart(v) < 0.0) {
29: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPart(v), (double)(-PetscImaginaryPart(v))));
30: } else {
31: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v)));
32: }
33: #else
34: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v));
35: #endif
36: }
37: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " constrained"));
38: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
39: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
40: } else {
41: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
42: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) {
43: PetscScalar v = array[i];
44: #if defined(PETSC_USE_COMPLEX)
45: if (PetscImaginaryPart(v) > 0.0) {
46: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v)));
47: } else if (PetscImaginaryPart(v) < 0.0) {
48: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPart(v), (double)(-PetscImaginaryPart(v))));
49: } else {
50: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v)));
51: }
52: #else
53: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v));
54: #endif
55: }
56: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
57: }
58: }
59: PetscCall(PetscViewerFlush(viewer));
60: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
61: PetscCall(VecRestoreArray(v, &array));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: /*@
66: PetscSectionVecView - View a vector, using the section to structure the values
68: Not Collective
70: Input Parameters:
71: + s - the organizing `PetscSection`
72: . v - the `Vec`
73: - viewer - the `PetscViewer`
75: Level: developer
77: .seealso: `PetscSection`, `PetscViewer`, `PetscSectionCreate()`, `VecSetValuesSection()`
78: @*/
79: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
80: {
81: PetscBool isascii;
82: PetscInt f;
84: PetscFunctionBegin;
87: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer));
89: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
90: if (isascii) {
91: const char *name;
93: PetscCall(PetscObjectGetName((PetscObject)v, &name));
94: if (s->numFields) {
95: PetscCall(PetscViewerASCIIPrintf(viewer, "%s with %" PetscInt_FMT " fields\n", name, s->numFields));
96: for (f = 0; f < s->numFields; ++f) {
97: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
98: PetscCall(PetscSectionVecView_ASCII(s->field[f], v, viewer));
99: }
100: } else {
101: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", name));
102: PetscCall(PetscSectionVecView_ASCII(s, v, viewer));
103: }
104: }
105: PetscFunctionReturn(PETSC_SUCCESS);
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;
130: PetscFunctionBegin;
133: PetscCall(VecGetArray(v, &baseArray));
134: *values = &baseArray[s->atlasOff[p]];
135: PetscCall(VecRestoreArray(v, &baseArray));
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: /*@
140: VecSetValuesSection - Sets all the values associated with a given point, according to the section, in the given `Vec`
142: Not Collective
144: Input Parameters:
145: + v - the `Vec`
146: . s - the organizing `PetscSection`
147: . point - the point
148: . values - the array of input values
149: - mode - the insertion mode, either `ADD_VALUES` or `INSERT_VALUES`
151: Level: developer
153: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecGetValuesSection()`
154: @*/
155: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, const PetscScalar values[], InsertMode mode)
156: {
157: PetscScalar *baseArray, *array;
158: const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
159: const PetscBool doInterior = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
160: const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
161: const PetscInt p = point - s->pStart;
162: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
163: PetscInt cDim = 0;
165: PetscFunctionBegin;
168: PetscCall(PetscSectionGetConstraintDof(s, point, &cDim));
169: PetscCall(VecGetArray(v, &baseArray));
170: array = &baseArray[s->atlasOff[p]];
171: if (!cDim && doInterior) {
172: if (orientation >= 0) {
173: const PetscInt dim = s->atlasDof[p];
174: PetscInt i;
176: if (doInsert) {
177: for (i = 0; i < dim; ++i) array[i] = values[i];
178: } else {
179: for (i = 0; i < dim; ++i) array[i] += values[i];
180: }
181: } else {
182: PetscInt offset = 0;
183: PetscInt j = -1, field, i;
185: for (field = 0; field < s->numFields; ++field) {
186: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
188: for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
189: offset += dim;
190: }
191: }
192: } else if (cDim) {
193: if (orientation >= 0) {
194: const PetscInt dim = s->atlasDof[p];
195: PetscInt cInd = 0, i;
196: const PetscInt *cDof;
198: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
199: if (doInsert) {
200: for (i = 0; i < dim; ++i) {
201: if ((cInd < cDim) && (i == cDof[cInd])) {
202: if (doBC) array[i] = values[i]; /* Constrained update */
203: ++cInd;
204: continue;
205: }
206: if (doInterior) array[i] = values[i]; /* Unconstrained update */
207: }
208: } else {
209: for (i = 0; i < dim; ++i) {
210: if ((cInd < cDim) && (i == cDof[cInd])) {
211: if (doBC) array[i] += values[i]; /* Constrained update */
212: ++cInd;
213: continue;
214: }
215: if (doInterior) array[i] += values[i]; /* Unconstrained update */
216: }
217: }
218: } else {
219: /* TODO This is broken for add and constrained update */
220: const PetscInt *cDof;
221: PetscInt offset = 0;
222: PetscInt cOffset = 0;
223: PetscInt j = 0, field;
225: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
226: for (field = 0; field < s->numFields; ++field) {
227: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
228: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
229: const PetscInt sDim = dim - tDim;
230: PetscInt cInd = 0, i, k;
232: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
233: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
234: ++cInd;
235: continue;
236: }
237: if (doInterior) array[j] = values[k]; /* Unconstrained update */
238: }
239: offset += dim;
240: cOffset += dim - tDim;
241: }
242: }
243: }
244: PetscCall(VecRestoreArray(v, &baseArray));
245: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
254: PetscCall(PetscSectionGetFieldComponents(section, field, &Nc));
255: for (p = pStart; p < pEnd; ++p) {
256: PetscInt gdof, fdof = 0;
258: PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
259: if (gdof > 0) PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
260: subSize += fdof;
261: }
262: PetscCall(PetscMalloc1(subSize, &subIndices));
263: for (p = pStart; p < pEnd; ++p) {
264: PetscInt gdof, goff;
266: PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
267: if (gdof > 0) {
268: PetscInt fdof, fc, f2, poff = 0;
270: PetscCall(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: PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof));
274: poff += fdof;
275: }
276: PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
277: for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
278: }
279: }
280: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)v), subSize, subIndices, PETSC_OWN_POINTER, is));
281: PetscCall(VecGetSubVector(v, *is, subv));
282: PetscCall(VecSetBlockSize(*subv, Nc));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
287: {
288: PetscFunctionBegin;
289: PetscCall(VecRestoreSubVector(v, *is, subv));
290: PetscCall(ISDestroy(is));
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /*@
295: PetscSectionVecNorm - Computes the vector norm of each field
297: Input Parameters:
298: + s - the local Section
299: . gs - the global section
300: . x - the vector
301: - type - one of `NORM_1`, `NORM_2`, `NORM_INFINITY`.
303: Output Parameter:
304: . val - the array of norms
306: Level: intermediate
308: .seealso: `VecNorm()`, `PetscSectionCreate()`
309: @*/
310: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
311: {
312: PetscInt Nf, f, pStart, pEnd;
314: PetscFunctionBegin;
318: PetscAssertPointer(val, 5);
319: PetscCall(PetscSectionGetNumFields(s, &Nf));
320: if (Nf < 2) PetscCall(VecNorm(x, type, val));
321: else {
322: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
323: for (f = 0; f < Nf; ++f) {
324: Vec subv;
325: IS is;
327: PetscCall(PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv));
328: PetscCall(VecNorm(subv, type, &val[f]));
329: PetscCall(PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv));
330: }
331: }
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }