Actual source code: vsection.c
petsc-3.5.4 2015-05-23
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: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
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: VecRestoreArray(v, &array);
66: return(0);
67: }
71: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
72: {
73: PetscBool isascii;
74: PetscInt f;
78: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)v), &viewer);}
81: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
82: if (isascii) {
83: const char *name;
85: PetscObjectGetName((PetscObject) v, &name);
86: if (s->numFields) {
87: PetscViewerASCIIPrintf(viewer, "%s with %d fields\n", name, s->numFields);
88: for (f = 0; f < s->numFields; ++f) {
89: PetscViewerASCIIPrintf(viewer, " field %d with %d components\n", f, s->numFieldComponents[f]);
90: PetscSectionVecView_ASCII(s->field[f], v, viewer);
91: }
92: } else {
93: PetscViewerASCIIPrintf(viewer, "%s\n", name);
94: PetscSectionVecView_ASCII(s, v, viewer);
95: }
96: }
97: return(0);
98: }
102: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
103: {
104: PetscScalar *baseArray;
105: const PetscInt p = point - s->pStart;
109: VecGetArray(v, &baseArray);
110: *values = &baseArray[s->atlasOff[p]];
111: VecRestoreArray(v, &baseArray);
112: return(0);
113: }
117: /*@C
118: VecSetValuesSection - Sets all the values associated with a given point, accoridng to the section, in the given Vec
120: Not collective
122: Input Parameters:
123: + v - the Vec
124: . s - the organizing PetscSection
125: . point - the point
126: . values - the array of input values
127: - mode - the insertion mode, either ADD_VALUES or INSERT_VALUES
129: Level: developer
131: Note: This is similar to MatSetValuesStencil(). The Fortran binding is
132: $
133: $ VecSetValuesSectionF90(vec, section, point, values, mode, ierr)
134: $
136: .seealso: PetscSection, PetscSectionCreate()
137: @*/
138: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
139: {
140: PetscScalar *baseArray, *array;
141: const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
142: const PetscBool doInterior = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
143: const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
144: const PetscInt p = point - s->pStart;
145: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
146: PetscInt cDim = 0;
147: PetscErrorCode ierr;
150: PetscSectionGetConstraintDof(s, point, &cDim);
151: VecGetArray(v, &baseArray);
152: array = &baseArray[s->atlasOff[p]];
153: if (!cDim && doInterior) {
154: if (orientation >= 0) {
155: const PetscInt dim = s->atlasDof[p];
156: PetscInt i;
158: if (doInsert) {
159: for (i = 0; i < dim; ++i) array[i] = values[i];
160: } else {
161: for (i = 0; i < dim; ++i) array[i] += values[i];
162: }
163: } else {
164: PetscInt offset = 0;
165: PetscInt j = -1, field, i;
167: for (field = 0; field < s->numFields; ++field) {
168: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
170: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
171: offset += dim;
172: }
173: }
174: } else if (cDim) {
175: if (orientation >= 0) {
176: const PetscInt dim = s->atlasDof[p];
177: PetscInt cInd = 0, i;
178: const PetscInt *cDof;
180: PetscSectionGetConstraintIndices(s, point, &cDof);
181: if (doInsert) {
182: for (i = 0; i < dim; ++i) {
183: if ((cInd < cDim) && (i == cDof[cInd])) {
184: if (doBC) array[i] = values[i]; /* Constrained update */
185: ++cInd;
186: continue;
187: }
188: if (doInterior) array[i] = values[i]; /* Unconstrained update */
189: }
190: } else {
191: for (i = 0; i < dim; ++i) {
192: if ((cInd < cDim) && (i == cDof[cInd])) {
193: if (doBC) array[i] += values[i]; /* Constrained update */
194: ++cInd;
195: continue;
196: }
197: if (doInterior) array[i] += values[i]; /* Unconstrained update */
198: }
199: }
200: } else {
201: /* TODO This is broken for add and constrained update */
202: const PetscInt *cDof;
203: PetscInt offset = 0;
204: PetscInt cOffset = 0;
205: PetscInt j = 0, field;
207: PetscSectionGetConstraintIndices(s, point, &cDof);
208: for (field = 0; field < s->numFields; ++field) {
209: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
210: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
211: const PetscInt sDim = dim - tDim;
212: PetscInt cInd = 0, i ,k;
214: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
215: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
216: if (doInterior) array[j] = values[k]; /* Unconstrained update */
217: }
218: offset += dim;
219: cOffset += dim - tDim;
220: }
221: }
222: }
223: VecRestoreArray(v, &baseArray);
224: return(0);
225: }