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: /*@
8: PetscSectionVecView - View a vector, using the section to structure the values
10: Collective
12: Input Parameters:
13: + s - the organizing `PetscSection`
14: . v - the `Vec`
15: - viewer - the `PetscViewer`
17: Level: developer
19: .seealso: `PetscSection`, `PetscViewer`, `PetscSectionCreate()`, `VecSetValuesSection()`, `PetscSectionArrayView()`
20: @*/
21: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
22: {
23: PetscBool isascii;
24: PetscScalar *array;
26: PetscFunctionBegin;
29: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer));
31: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
32: if (isascii) {
33: const char *name;
35: PetscCall(PetscObjectGetName((PetscObject)v, &name));
36: if (s->numFields) {
37: PetscCall(PetscViewerASCIIPrintf(viewer, "%s with %" PetscInt_FMT " fields\n", name, s->numFields));
38: for (PetscInt f = 0; f < s->numFields; ++f) {
39: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
40: PetscCall(VecGetArray(v, &array));
41: PetscCall(PetscSectionArrayView_ASCII_Internal(s->field[f], array, PETSC_SCALAR, viewer));
42: PetscCall(VecRestoreArray(v, &array));
43: }
44: } else {
45: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", name));
46: PetscCall(VecGetArray(v, &array));
47: PetscCall(PetscSectionArrayView_ASCII_Internal(s, array, PETSC_SCALAR, viewer));
48: PetscCall(VecRestoreArray(v, &array));
49: }
50: }
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*@C
55: VecGetValuesSection - Gets all the values associated with a given point, according to the section, in the given `Vec`
57: Not Collective
59: Input Parameters:
60: + v - the `Vec`
61: . s - the organizing `PetscSection`
62: - point - the point
64: Output Parameter:
65: . values - the array of output values
67: Level: developer
69: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecSetValuesSection()`
70: @*/
71: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar *values[])
72: {
73: PetscScalar *baseArray;
74: const PetscInt p = point - s->pStart;
76: PetscFunctionBegin;
79: PetscCall(VecGetArray(v, &baseArray));
80: *values = &baseArray[s->atlasOff[p]];
81: PetscCall(VecRestoreArray(v, &baseArray));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@
86: VecSetValuesSection - Sets all the values associated with a given point, according to the section, in the given `Vec`
88: Not Collective
90: Input Parameters:
91: + v - the `Vec`
92: . s - the organizing `PetscSection`
93: . point - the point
94: . values - the array of input values
95: - mode - the insertion mode, either `ADD_VALUES` or `INSERT_VALUES`
97: Level: developer
99: .seealso: `PetscSection`, `PetscSectionCreate()`, `VecGetValuesSection()`
100: @*/
101: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, const PetscScalar values[], InsertMode mode)
102: {
103: PetscScalar *baseArray, *array;
104: const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
105: const PetscBool doInterior = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
106: const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
107: const PetscInt p = point - s->pStart;
108: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
109: PetscInt cDim = 0;
111: PetscFunctionBegin;
114: PetscCall(PetscSectionGetConstraintDof(s, point, &cDim));
115: PetscCall(VecGetArray(v, &baseArray));
116: array = &baseArray[s->atlasOff[p]];
117: if (!cDim && doInterior) {
118: if (orientation >= 0) {
119: const PetscInt dim = s->atlasDof[p];
120: PetscInt i;
122: if (doInsert) {
123: for (i = 0; i < dim; ++i) array[i] = values[i];
124: } else {
125: for (i = 0; i < dim; ++i) array[i] += values[i];
126: }
127: } else {
128: PetscInt offset = 0;
129: PetscInt j = -1, field, i;
131: for (field = 0; field < s->numFields; ++field) {
132: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
134: for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
135: offset += dim;
136: }
137: }
138: } else if (cDim) {
139: if (orientation >= 0) {
140: const PetscInt dim = s->atlasDof[p];
141: PetscInt cInd = 0, i;
142: const PetscInt *cDof;
144: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
145: if (doInsert) {
146: for (i = 0; i < dim; ++i) {
147: if ((cInd < cDim) && (i == cDof[cInd])) {
148: if (doBC) array[i] = values[i]; /* Constrained update */
149: ++cInd;
150: continue;
151: }
152: if (doInterior) array[i] = values[i]; /* Unconstrained update */
153: }
154: } else {
155: for (i = 0; i < dim; ++i) {
156: if ((cInd < cDim) && (i == cDof[cInd])) {
157: if (doBC) array[i] += values[i]; /* Constrained update */
158: ++cInd;
159: continue;
160: }
161: if (doInterior) array[i] += values[i]; /* Unconstrained update */
162: }
163: }
164: } else {
165: /* TODO This is broken for add and constrained update */
166: const PetscInt *cDof;
167: PetscInt offset = 0;
168: PetscInt cOffset = 0;
169: PetscInt j = 0, field;
171: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
172: for (field = 0; field < s->numFields; ++field) {
173: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
174: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
175: const PetscInt sDim = dim - tDim;
176: PetscInt cInd = 0, i, k;
178: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
179: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
180: ++cInd;
181: continue;
182: }
183: if (doInterior) array[j] = values[k]; /* Unconstrained update */
184: }
185: offset += dim;
186: cOffset += dim - tDim;
187: }
188: }
189: }
190: PetscCall(VecRestoreArray(v, &baseArray));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
195: {
196: PetscInt *subIndices;
197: PetscInt Nc, subSize = 0, subOff = 0;
199: PetscFunctionBegin;
200: PetscCall(PetscSectionGetFieldComponents(section, field, &Nc));
201: for (PetscInt p = pStart; p < pEnd; ++p) {
202: PetscInt gdof, fdof = 0;
204: PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
205: if (gdof > 0) PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
206: subSize += fdof;
207: }
208: PetscCall(PetscMalloc1(subSize, &subIndices));
209: for (PetscInt p = pStart; p < pEnd; ++p) {
210: PetscInt gdof, goff;
212: PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
213: if (gdof > 0) {
214: PetscInt fdof, fc, f2, poff = 0;
216: PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
217: /* Can get rid of this loop by storing field information in the global section */
218: for (f2 = 0; f2 < field; ++f2) {
219: PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof));
220: poff += fdof;
221: }
222: PetscCall(PetscSectionGetFieldDof(section, p, field, &fdof));
223: for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
224: }
225: }
226: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)v), subSize, subIndices, PETSC_OWN_POINTER, is));
227: PetscCall(VecGetSubVector(v, *is, subv));
228: PetscCall(VecSetBlockSize(*subv, Nc));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
233: {
234: PetscFunctionBegin;
235: PetscCall(VecRestoreSubVector(v, *is, subv));
236: PetscCall(ISDestroy(is));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@
241: PetscSectionVecNorm - Computes the vector norm of each field
243: Input Parameters:
244: + s - the local Section
245: . gs - the global section
246: . x - the vector
247: - type - one of `NORM_1`, `NORM_2`, `NORM_INFINITY`.
249: Output Parameter:
250: . val - the array of norms
252: Level: intermediate
254: .seealso: `VecNorm()`, `PetscSectionCreate()`
255: @*/
256: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
257: {
258: PetscInt Nf, f, pStart, pEnd;
260: PetscFunctionBegin;
264: PetscAssertPointer(val, 5);
265: PetscCall(PetscSectionGetNumFields(s, &Nf));
266: if (Nf < 2) PetscCall(VecNorm(x, type, val));
267: else {
268: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
269: for (f = 0; f < Nf; ++f) {
270: Vec subv;
271: IS is;
273: PetscCall(PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv));
274: PetscCall(VecNorm(subv, type, &val[f]));
275: PetscCall(PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv));
276: }
277: }
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }