Actual source code: vsection.c
petsc-3.7.7 2017-09-25
1: /*
2: This file contains routines for section object operations on Vecs
3: */
4: #include <petsc/private/isimpl.h> /*I "petscvec.h" I*/
5: #include <petsc/private/vecimpl.h> /*I "petscvec.h" I*/
9: PetscErrorCode PetscSectionVecView_ASCII(PetscSection s, Vec v, PetscViewer viewer)
10: {
11: PetscScalar *array;
12: PetscInt p, i;
13: PetscMPIInt rank;
17: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
18: VecGetArray(v, &array);
19: PetscViewerASCIIPushSynchronized(viewer);
20: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
21: for (p = 0; p < s->pEnd - s->pStart; ++p) {
22: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
23: PetscInt b;
25: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
26: for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
27: PetscScalar v = array[i];
28: #if defined(PETSC_USE_COMPLEX)
29: if (PetscImaginaryPart(v) > 0.0) {
30: PetscViewerASCIISynchronizedPrintf(viewer," %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));
31: } else if (PetscImaginaryPart(v) < 0.0) {
32: PetscViewerASCIISynchronizedPrintf(viewer," %g - %g i", (double)PetscRealPart(v),(double)(-PetscImaginaryPart(v)));
33: } else {
34: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));
35: }
36: #else
37: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);
38: #endif
39: }
40: PetscViewerASCIISynchronizedPrintf(viewer, " constrained");
41: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
42: PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
43: }
44: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
45: } else {
46: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
47: for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
48: PetscScalar v = array[i];
49: #if defined(PETSC_USE_COMPLEX)
50: if (PetscImaginaryPart(v) > 0.0) {
51: PetscViewerASCIISynchronizedPrintf(viewer," %g + %g i", (double)PetscRealPart(v), (double)PetscImaginaryPart(v));
52: } else if (PetscImaginaryPart(v) < 0.0) {
53: PetscViewerASCIISynchronizedPrintf(viewer," %g - %g i", (double)PetscRealPart(v),(double)(-PetscImaginaryPart(v)));
54: } else {
55: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPart(v));
56: }
57: #else
58: PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)v);
59: #endif
60: }
61: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
62: }
63: }
64: PetscViewerFlush(viewer);
65: PetscViewerASCIIPopSynchronized(viewer);
66: VecRestoreArray(v, &array);
67: return(0);
68: }
72: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
73: {
74: PetscBool isascii;
75: PetscInt f;
79: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);}
82: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
83: if (isascii) {
84: const char *name;
86: PetscObjectGetName((PetscObject) v, &name);
87: if (s->numFields) {
88: PetscViewerASCIIPrintf(viewer, "%s with %D fields\n", name, s->numFields);
89: for (f = 0; f < s->numFields; ++f) {
90: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
91: PetscSectionVecView_ASCII(s->field[f], v, viewer);
92: }
93: } else {
94: PetscViewerASCIIPrintf(viewer, "%s\n", name);
95: PetscSectionVecView_ASCII(s, v, viewer);
96: }
97: }
98: return(0);
99: }
103: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
104: {
105: PetscScalar *baseArray;
106: const PetscInt p = point - s->pStart;
110: VecGetArray(v, &baseArray);
111: *values = &baseArray[s->atlasOff[p]];
112: VecRestoreArray(v, &baseArray);
113: return(0);
114: }
118: /*@C
119: VecSetValuesSection - Sets all the values associated with a given point, according to the section, in the given Vec
121: Not collective
123: Input Parameters:
124: + v - the Vec
125: . s - the organizing PetscSection
126: . point - the point
127: . values - the array of input values
128: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES
130: Level: developer
132: Note: This is similar to MatSetValuesStencil(). The Fortran binding is
133: $
134: $ VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
135: $
137: .seealso: PetscSection, PetscSectionCreate()
138: @*/
139: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
140: {
141: PetscScalar *baseArray, *array;
142: const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
143: const PetscBool doInterior = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
144: const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
145: const PetscInt p = point - s->pStart;
146: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
147: PetscInt cDim = 0;
148: PetscErrorCode ierr;
151: PetscSectionGetConstraintDof(s, point, &cDim);
152: VecGetArray(v, &baseArray);
153: array = &baseArray[s->atlasOff[p]];
154: if (!cDim && doInterior) {
155: if (orientation >= 0) {
156: const PetscInt dim = s->atlasDof[p];
157: PetscInt i;
159: if (doInsert) {
160: for (i = 0; i < dim; ++i) array[i] = values[i];
161: } else {
162: for (i = 0; i < dim; ++i) array[i] += values[i];
163: }
164: } else {
165: PetscInt offset = 0;
166: PetscInt j = -1, field, i;
168: for (field = 0; field < s->numFields; ++field) {
169: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
171: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
172: offset += dim;
173: }
174: }
175: } else if (cDim) {
176: if (orientation >= 0) {
177: const PetscInt dim = s->atlasDof[p];
178: PetscInt cInd = 0, i;
179: const PetscInt *cDof;
181: PetscSectionGetConstraintIndices(s, point, &cDof);
182: if (doInsert) {
183: for (i = 0; i < dim; ++i) {
184: if ((cInd < cDim) && (i == cDof[cInd])) {
185: if (doBC) array[i] = values[i]; /* Constrained update */
186: ++cInd;
187: continue;
188: }
189: if (doInterior) array[i] = values[i]; /* Unconstrained update */
190: }
191: } else {
192: for (i = 0; i < dim; ++i) {
193: if ((cInd < cDim) && (i == cDof[cInd])) {
194: if (doBC) array[i] += values[i]; /* Constrained update */
195: ++cInd;
196: continue;
197: }
198: if (doInterior) array[i] += values[i]; /* Unconstrained update */
199: }
200: }
201: } else {
202: /* TODO This is broken for add and constrained update */
203: const PetscInt *cDof;
204: PetscInt offset = 0;
205: PetscInt cOffset = 0;
206: PetscInt j = 0, field;
208: PetscSectionGetConstraintIndices(s, point, &cDof);
209: for (field = 0; field < s->numFields; ++field) {
210: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
211: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
212: const PetscInt sDim = dim - tDim;
213: PetscInt cInd = 0, i ,k;
215: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
216: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
217: if (doInterior) array[j] = values[k]; /* Unconstrained update */
218: }
219: offset += dim;
220: cOffset += dim - tDim;
221: }
222: }
223: }
224: VecRestoreArray(v, &baseArray);
225: return(0);
226: }
230: PetscErrorCode PetscSectionGetField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
231: {
232: PetscInt *subIndices;
233: PetscInt Nc, subSize = 0, subOff = 0, p;
237: PetscSectionGetFieldComponents(section, field, &Nc);
238: for (p = pStart; p < pEnd; ++p) {
239: PetscInt gdof, fdof = 0;
241: PetscSectionGetDof(sectionGlobal, p, &gdof);
242: if (gdof > 0) {PetscSectionGetFieldDof(section, p, field, &fdof);}
243: subSize += fdof;
244: }
245: PetscMalloc1(subSize, &subIndices);
246: for (p = pStart; p < pEnd; ++p) {
247: PetscInt gdof, goff;
249: PetscSectionGetDof(sectionGlobal, p, &gdof);
250: if (gdof > 0) {
251: PetscInt fdof, fc, f2, poff = 0;
253: PetscSectionGetOffset(sectionGlobal, p, &goff);
254: /* Can get rid of this loop by storing field information in the global section */
255: for (f2 = 0; f2 < field; ++f2) {
256: PetscSectionGetFieldDof(section, p, f2, &fdof);
257: poff += fdof;
258: }
259: PetscSectionGetFieldDof(section, p, field, &fdof);
260: for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff+poff+fc;
261: }
262: }
263: ISCreateGeneral(PetscObjectComm((PetscObject) v), subSize, subIndices, PETSC_OWN_POINTER, is);
264: VecGetSubVector(v, *is, subv);
265: VecSetBlockSize(*subv, Nc);
266: return(0);
267: }
271: PetscErrorCode PetscSectionRestoreField_Internal(PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
272: {
276: VecRestoreSubVector(v, *is, subv);
277: ISDestroy(is);
278: return(0);
279: }
283: /*@C
284: PetscSectionVecNorm - Computes the vector norm, separated into field components.
286: Input Parameters:
287: + s - the local Section
288: . gs - the global section
289: . x - the vector
290: - type - one of NORM_1, NORM_2, NORM_INFINITY.
292: Output Parameter:
293: . val - the array of norms
295: Level: intermediate
297: .seealso: VecNorm(), PetscSectionCreate()
298: @*/
299: PetscErrorCode PetscSectionVecNorm(PetscSection s, PetscSection gs, Vec x, NormType type, PetscReal val[])
300: {
301: PetscInt Nf, f, pStart, pEnd;
309: PetscSectionGetNumFields(s, &Nf);
310: if (Nf < 2) {VecNorm(x, type, val);}
311: else {
312: PetscSectionGetChart(s, &pStart, &pEnd);
313: for (f = 0; f < Nf; ++f) {
314: Vec subv;
315: IS is;
317: PetscSectionGetField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
318: VecNorm(subv, type, &val[f]);
319: PetscSectionRestoreField_Internal(s, gs, x, f, pStart, pEnd, &is, &subv);
320: }
321: }
322: return(0);
323: }