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 ") dim %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 ") dim %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: /*@C
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: Fortran Notes:
154: This is similar to `MatSetValuesStencil()`. The binding is
155: $ VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
157: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecGetValuesSection()`
158: @*/
159: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, const 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;
169: PetscFunctionBegin;
172: PetscCall(PetscSectionGetConstraintDof(s, point, &cDim));
173: PetscCall(VecGetArray(v, &baseArray));
174: array = &baseArray[s->atlasOff[p]];
175: if (!cDim && doInterior) {
176: if (orientation >= 0) {
177: const PetscInt dim = s->atlasDof[p];
178: PetscInt i;
180: if (doInsert) {
181: for (i = 0; i < dim; ++i) array[i] = values[i];
182: } else {
183: for (i = 0; i < dim; ++i) array[i] += values[i];
184: }
185: } else {
186: PetscInt offset = 0;
187: PetscInt j = -1, field, i;
189: for (field = 0; field < s->numFields; ++field) {
190: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
192: for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
193: offset += dim;
194: }
195: }
196: } else if (cDim) {
197: if (orientation >= 0) {
198: const PetscInt dim = s->atlasDof[p];
199: PetscInt cInd = 0, i;
200: const PetscInt *cDof;
202: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
203: if (doInsert) {
204: for (i = 0; i < dim; ++i) {
205: if ((cInd < cDim) && (i == cDof[cInd])) {
206: if (doBC) array[i] = values[i]; /* Constrained update */
207: ++cInd;
208: continue;
209: }
210: if (doInterior) array[i] = values[i]; /* Unconstrained update */
211: }
212: } else {
213: for (i = 0; i < dim; ++i) {
214: if ((cInd < cDim) && (i == cDof[cInd])) {
215: if (doBC) array[i] += values[i]; /* Constrained update */
216: ++cInd;
217: continue;
218: }
219: if (doInterior) array[i] += values[i]; /* Unconstrained update */
220: }
221: }
222: } else {
223: /* TODO This is broken for add and constrained update */
224: const PetscInt *cDof;
225: PetscInt offset = 0;
226: PetscInt cOffset = 0;
227: PetscInt j = 0, field;
229: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
230: for (field = 0; field < s->numFields; ++field) {
231: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
232: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
233: const PetscInt sDim = dim - tDim;
234: PetscInt cInd = 0, i, k;
236: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
237: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
238: ++cInd;
239: continue;
240: }
241: if (doInterior) array[j] = values[k]; /* Unconstrained update */
242: }
243: offset += dim;
244: cOffset += dim - tDim;
245: }
246: }
247: }
248: PetscCall(VecRestoreArray(v, &baseArray));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
253: {
254: PetscInt *subIndices;
255: PetscInt Nc, subSize = 0, subOff = 0, p;
257: PetscFunctionBegin;
258: PetscCall(PetscSectionGetFieldComponents(section, field, &Nc));
259: for (p = pStart; p < pEnd; ++p) {
260: PetscInt gdof, fdof = 0;
262: PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
263: if (gdof > 0) PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
264: subSize += fdof;
265: }
266: PetscCall(PetscMalloc1(subSize, &subIndices));
267: for (p = pStart; p < pEnd; ++p) {
268: PetscInt gdof, goff;
270: PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
271: if (gdof > 0) {
272: PetscInt fdof, fc, f2, poff = 0;
274: PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
275: /* Can get rid of this loop by storing field information in the global section */
276: for (f2 = 0; f2 < field; ++f2) {
277: PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof));
278: poff += fdof;
279: }
280: PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
281: for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
282: }
283: }
284: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)v), subSize, subIndices, PETSC_OWN_POINTER, is));
285: PetscCall(VecGetSubVector(v, *is, subv));
286: PetscCall(VecSetBlockSize(*subv, Nc));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
291: {
292: PetscFunctionBegin;
293: PetscCall(VecRestoreSubVector(v, *is, subv));
294: PetscCall(ISDestroy(is));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: /*@C
299: PetscSectionVecNorm - Computes the vector norm, separated into field components.
301: Input Parameters:
302: + s - the local Section
303: . gs - the global section
304: . x - the vector
305: - type - one of `NORM_1`, `NORM_2`, `NORM_INFINITY`.
307: Output Parameter:
308: . val - the array of norms
310: Level: intermediate
312: .seealso: `VecNorm()`, `PetscSectionCreate()`
313: @*/
314: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
315: {
316: PetscInt Nf, f, pStart, pEnd;
318: PetscFunctionBegin;
322: PetscAssertPointer(val, 5);
323: PetscCall(PetscSectionGetNumFields(s, &Nf));
324: if (Nf < 2) PetscCall(VecNorm(x, type, val));
325: else {
326: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
327: for (f = 0; f < Nf; ++f) {
328: Vec subv;
329: IS is;
331: PetscCall(PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv));
332: PetscCall(VecNorm(subv, type, &val[f]));
333: PetscCall(PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv));
334: }
335: }
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }