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