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
PetscSectionVecView_ASCII(PetscSection s,Vec v,PetscViewer viewer)7 static PetscErrorCode PetscSectionVecView_ASCII(PetscSection s, Vec v, PetscViewer viewer)
8 {
9 PetscScalar *array;
10 PetscInt p, i;
11 PetscMPIInt rank;
12 PetscErrorCode ierr;
13
14 PetscFunctionBegin;
15 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);CHKERRQ(ierr);
16 ierr = VecGetArray(v, &array);CHKERRQ(ierr);
17 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
18 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);CHKERRQ(ierr);
19 for (p = 0; p < s->pEnd - s->pStart; ++p) {
20 if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
21 PetscInt b;
22
23 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);CHKERRQ(ierr);
24 for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
25 PetscScalar v = array[i];
26 #if defined(PETSC_USE_COMPLEX)
27 if (PetscImaginaryPart(v) > 0.0) {
28 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));CHKERRQ(ierr);
29 } else if (PetscImaginaryPart(v) < 0.0) {
30 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %g - %g i", (double)PetscRealPart(v),(double)(-PetscImaginaryPart(v)));CHKERRQ(ierr);
31 } else {
32 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));CHKERRQ(ierr);
33 }
34 #else
35 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);CHKERRQ(ierr);
36 #endif
37 }
38 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " constrained");CHKERRQ(ierr);
39 for (b = 0; b < s->bc->atlasDof[p]; ++b) {
40 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);CHKERRQ(ierr);
41 }
42 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "\n");CHKERRQ(ierr);
43 } else {
44 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);CHKERRQ(ierr);
45 for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
46 PetscScalar v = array[i];
47 #if defined(PETSC_USE_COMPLEX)
48 if (PetscImaginaryPart(v) > 0.0) {
49 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));CHKERRQ(ierr);
50 } else if (PetscImaginaryPart(v) < 0.0) {
51 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %g - %g i", (double)PetscRealPart(v),(double)(-PetscImaginaryPart(v)));CHKERRQ(ierr);
52 } else {
53 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));CHKERRQ(ierr);
54 }
55 #else
56 ierr = PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);CHKERRQ(ierr);
57 #endif
58 }
59 ierr = PetscViewerASCIISynchronizedPrintf(viewer, "\n");CHKERRQ(ierr);
60 }
61 }
62 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
63 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
64 ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
65 PetscFunctionReturn(0);
66 }
67
68 /*@
69 PetscSectionVecView - View a vector, using the section to structure the values
70
71 Not collective
72
73 Input Parameters:
74 + s - the organizing PetscSection
75 . v - the Vec
76 - viewer - the PetscViewer
77
78 Level: developer
79
80 .seealso: PetscSection, PetscSectionCreate(), VecSetValuesSection()
81 @*/
PetscSectionVecView(PetscSection s,Vec v,PetscViewer viewer)82 PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
83 {
84 PetscBool isascii;
85 PetscInt f;
86 PetscErrorCode ierr;
87
88 PetscFunctionBegin;
89 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
90 PetscValidHeaderSpecific(v, VEC_CLASSID, 2);
91 if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);CHKERRQ(ierr);}
92 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 3);
93 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
94 if (isascii) {
95 const char *name;
96
97 ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr);
98 if (s->numFields) {
99 ierr = PetscViewerASCIIPrintf(viewer, "%s with %D fields\n", name, s->numFields);CHKERRQ(ierr);
100 for (f = 0; f < s->numFields; ++f) {
101 ierr = PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);CHKERRQ(ierr);
102 ierr = PetscSectionVecView_ASCII(s->field[f], v, viewer);CHKERRQ(ierr);
103 }
104 } else {
105 ierr = PetscViewerASCIIPrintf(viewer, "%s\n", name);CHKERRQ(ierr);
106 ierr = PetscSectionVecView_ASCII(s, v, viewer);CHKERRQ(ierr);
107 }
108 }
109 PetscFunctionReturn(0);
110 }
111
112 /*@C
113 VecGetValuesSection - Gets all the values associated with a given point, according to the section, in the given Vec
114
115 Not collective
116
117 Input Parameters:
118 + v - the Vec
119 . s - the organizing PetscSection
120 - point - the point
121
122 Output Parameter:
123 . values - the array of output values
124
125 Level: developer
126
127 .seealso: PetscSection, PetscSectionCreate(), VecSetValuesSection()
128 @*/
VecGetValuesSection(Vec v,PetscSection s,PetscInt point,PetscScalar ** values)129 PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
130 {
131 PetscScalar *baseArray;
132 const PetscInt p = point - s->pStart;
133 PetscErrorCode ierr;
134
135 PetscFunctionBegin;
136 PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
137 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2);
138 ierr = VecGetArray(v, &baseArray);CHKERRQ(ierr);
139 *values = &baseArray[s->atlasOff[p]];
140 ierr = VecRestoreArray(v, &baseArray);CHKERRQ(ierr);
141 PetscFunctionReturn(0);
142 }
143
144 /*@C
145 VecSetValuesSection - Sets all the values associated with a given point, according to the section, in the given Vec
146
147 Not collective
148
149 Input Parameters:
150 + v - the Vec
151 . s - the organizing PetscSection
152 . point - the point
153 . values - the array of input values
154 - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES
155
156 Level: developer
157
158 Note: This is similar to MatSetValuesStencil(). The Fortran binding is
159 $
160 $ VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
161 $
162
163 .seealso: PetscSection, PetscSectionCreate(), VecGetValuesSection()
164 @*/
VecSetValuesSection(Vec v,PetscSection s,PetscInt point,PetscScalar values[],InsertMode mode)165 PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
166 {
167 PetscScalar *baseArray, *array;
168 const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
169 const PetscBool doInterior = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
170 const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
171 const PetscInt p = point - s->pStart;
172 const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
173 PetscInt cDim = 0;
174 PetscErrorCode ierr;
175
176 PetscFunctionBegin;
177 PetscValidHeaderSpecific(v, VEC_CLASSID, 1);
178 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 2);
179 ierr = PetscSectionGetConstraintDof(s, point, &cDim);CHKERRQ(ierr);
180 ierr = VecGetArray(v, &baseArray);CHKERRQ(ierr);
181 array = &baseArray[s->atlasOff[p]];
182 if (!cDim && doInterior) {
183 if (orientation >= 0) {
184 const PetscInt dim = s->atlasDof[p];
185 PetscInt i;
186
187 if (doInsert) {
188 for (i = 0; i < dim; ++i) array[i] = values[i];
189 } else {
190 for (i = 0; i < dim; ++i) array[i] += values[i];
191 }
192 } else {
193 PetscInt offset = 0;
194 PetscInt j = -1, field, i;
195
196 for (field = 0; field < s->numFields; ++field) {
197 const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
198
199 for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
200 offset += dim;
201 }
202 }
203 } else if (cDim) {
204 if (orientation >= 0) {
205 const PetscInt dim = s->atlasDof[p];
206 PetscInt cInd = 0, i;
207 const PetscInt *cDof;
208
209 ierr = PetscSectionGetConstraintIndices(s, point, &cDof);CHKERRQ(ierr);
210 if (doInsert) {
211 for (i = 0; i < dim; ++i) {
212 if ((cInd < cDim) && (i == cDof[cInd])) {
213 if (doBC) array[i] = values[i]; /* Constrained update */
214 ++cInd;
215 continue;
216 }
217 if (doInterior) array[i] = values[i]; /* Unconstrained update */
218 }
219 } else {
220 for (i = 0; i < dim; ++i) {
221 if ((cInd < cDim) && (i == cDof[cInd])) {
222 if (doBC) array[i] += values[i]; /* Constrained update */
223 ++cInd;
224 continue;
225 }
226 if (doInterior) array[i] += values[i]; /* Unconstrained update */
227 }
228 }
229 } else {
230 /* TODO This is broken for add and constrained update */
231 const PetscInt *cDof;
232 PetscInt offset = 0;
233 PetscInt cOffset = 0;
234 PetscInt j = 0, field;
235
236 ierr = PetscSectionGetConstraintIndices(s, point, &cDof);CHKERRQ(ierr);
237 for (field = 0; field < s->numFields; ++field) {
238 const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
239 const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
240 const PetscInt sDim = dim - tDim;
241 PetscInt cInd = 0, i ,k;
242
243 for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
244 if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
245 if (doInterior) array[j] = values[k]; /* Unconstrained update */
246 }
247 offset += dim;
248 cOffset += dim - tDim;
249 }
250 }
251 }
252 ierr = VecRestoreArray(v, &baseArray);CHKERRQ(ierr);
253 PetscFunctionReturn(0);
254 }
255
PetscSectionGetField_Internal(PetscSection section,PetscSection sectionGlobal,Vec v,PetscInt field,PetscInt pStart,PetscInt pEnd,IS * is,Vec * subv)256 PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
257 {
258 PetscInt *subIndices;
259 PetscInt Nc, subSize = 0, subOff = 0, p;
260 PetscErrorCode ierr;
261
262 PetscFunctionBegin;
263 ierr = PetscSectionGetFieldComponents(section, field, &Nc);CHKERRQ(ierr);
264 for (p = pStart; p < pEnd; ++p) {
265 PetscInt gdof, fdof = 0;
266
267 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
268 if (gdof > 0) {ierr = PetscSectionGetFieldDof(section, p, field, &fdof);CHKERRQ(ierr);}
269 subSize += fdof;
270 }
271 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr);
272 for (p = pStart; p < pEnd; ++p) {
273 PetscInt gdof, goff;
274
275 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
276 if (gdof > 0) {
277 PetscInt fdof, fc, f2, poff = 0;
278
279 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
280 /* Can get rid of this loop by storing field information in the global section */
281 for (f2 = 0; f2 < field; ++f2) {
282 ierr = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr);
283 poff += fdof;
284 }
285 ierr = PetscSectionGetFieldDof(section, p, field, &fdof);CHKERRQ(ierr);
286 for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff+poff+fc;
287 }
288 }
289 ierr = ISCreateGeneral(PetscObjectComm((PetscObject) v), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr);
290 ierr = VecGetSubVector(v, *is, subv);CHKERRQ(ierr);
291 ierr = VecSetBlockSize(*subv, Nc);CHKERRQ(ierr);
292 PetscFunctionReturn(0);
293 }
294
PetscSectionRestoreField_Internal(PetscSection section,PetscSection sectionGlobal,Vec v,PetscInt field,PetscInt pStart,PetscInt pEnd,IS * is,Vec * subv)295 PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
296 {
297 PetscErrorCode ierr;
298
299 PetscFunctionBegin;
300 ierr = VecRestoreSubVector(v, *is, subv);CHKERRQ(ierr);
301 ierr = ISDestroy(is);CHKERRQ(ierr);
302 PetscFunctionReturn(0);
303 }
304
305 /*@C
306 PetscSectionVecNorm - Computes the vector norm, separated into field components.
307
308 Input Parameters:
309 + s - the local Section
310 . gs - the global section
311 . x - the vector
312 - type - one of NORM_1, NORM_2, NORM_INFINITY.
313
314 Output Parameter:
315 . val - the array of norms
316
317 Level: intermediate
318
319 .seealso: VecNorm(), PetscSectionCreate()
320 @*/
PetscSectionVecNorm(PetscSection s,PetscSection gs,Vec x,NormType type,PetscReal val[])321 PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
322 {
323 PetscInt Nf, f, pStart, pEnd;
324 PetscErrorCode ierr;
325
326 PetscFunctionBegin;
327 PetscValidHeaderSpecific(s, PETSC_SECTION_CLASSID, 1);
328 PetscValidHeaderSpecific(gs, PETSC_SECTION_CLASSID, 2);
329 PetscValidHeaderSpecific(x, VEC_CLASSID, 3);
330 PetscValidPointer(val, 5);
331 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr);
332 if (Nf < 2) {ierr = VecNorm(x, type, val);CHKERRQ(ierr);}
333 else {
334 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
335 for (f = 0; f < Nf; ++f) {
336 Vec subv;
337 IS is;
338
339 ierr = PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);CHKERRQ(ierr);
340 ierr = VecNorm(subv, type, &val[f]);CHKERRQ(ierr);
341 ierr = PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);CHKERRQ(ierr);
342 }
343 }
344 PetscFunctionReturn(0);
345 }
346