Actual source code: vsection.c
petsc-3.4.5 2014-06-29
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: if (s->atlasLayout.numDof != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle %d dof in a uniform section", s->atlasLayout.numDof);
18: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
19: VecGetArray(v, &array);
20: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
21: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
22: for (p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
23: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
24: PetscInt b;
26: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d", p+s->atlasLayout.pStart, s->atlasDof[p], s->atlasOff[p]);
27: for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
28: PetscScalar v = array[i];
29: #if defined(PETSC_USE_COMPLEX)
30: if (PetscImaginaryPart(v) > 0.0) {
31: PetscViewerASCIISynchronizedPrintf(viewer," %G + %G i", PetscRealPart(v), PetscImaginaryPart(v));
32: } else if (PetscImaginaryPart(v) < 0.0) {
33: PetscViewerASCIISynchronizedPrintf(viewer," %G - %G i", PetscRealPart(v),-PetscImaginaryPart(v));
34: } else {
35: PetscViewerASCIISynchronizedPrintf(viewer, " %G", PetscRealPart(v));
36: }
37: #else
38: PetscViewerASCIISynchronizedPrintf(viewer, " %G", v);
39: #endif
40: }
41: PetscViewerASCIISynchronizedPrintf(viewer, " constrained");
42: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
43: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
44: }
45: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
46: } else {
47: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d", p+s->atlasLayout.pStart, s->atlasDof[p], s->atlasOff[p]);
48: for (i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
49: PetscScalar v = array[i];
50: #if defined(PETSC_USE_COMPLEX)
51: if (PetscImaginaryPart(v) > 0.0) {
52: PetscViewerASCIISynchronizedPrintf(viewer," %G + %G i", PetscRealPart(v), PetscImaginaryPart(v));
53: } else if (PetscImaginaryPart(v) < 0.0) {
54: PetscViewerASCIISynchronizedPrintf(viewer," %G - %G i", PetscRealPart(v),-PetscImaginaryPart(v));
55: } else {
56: PetscViewerASCIISynchronizedPrintf(viewer, " %G", PetscRealPart(v));
57: }
58: #else
59: PetscViewerASCIISynchronizedPrintf(viewer, " %G", v);
60: #endif
61: }
62: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
63: }
64: }
65: PetscViewerFlush(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->atlasLayout.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, accoridng 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()
134: .seealso: PetscSection, PetscSectionCreate()
135: @*/
136: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
137: {
138: PetscScalar *baseArray, *array;
139: const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
140: const PetscBool doInterior = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_VALUES || mode == ADD_VALUES ? PETSC_TRUE : PETSC_FALSE;
141: const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES || mode == INSERT_BC_VALUES || mode == ADD_BC_VALUES ? PETSC_TRUE : PETSC_FALSE;
142: const PetscInt p = point - s->atlasLayout.pStart;
143: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
144: PetscInt cDim = 0;
145: PetscErrorCode ierr;
148: PetscSectionGetConstraintDof(s, p, &cDim);
149: VecGetArray(v, &baseArray);
150: array = &baseArray[s->atlasOff[p]];
151: if (!cDim && doInterior) {
152: if (orientation >= 0) {
153: const PetscInt dim = s->atlasDof[p];
154: PetscInt i;
156: if (doInsert) {
157: for (i = 0; i < dim; ++i) array[i] = values[i];
158: } else {
159: for (i = 0; i < dim; ++i) array[i] += values[i];
160: }
161: } else {
162: PetscInt offset = 0;
163: PetscInt j = -1, field, i;
165: for (field = 0; field < s->numFields; ++field) {
166: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
168: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
169: offset += dim;
170: }
171: }
172: } else if (cDim) {
173: if (orientation >= 0) {
174: const PetscInt dim = s->atlasDof[p];
175: PetscInt cInd = 0, i;
176: const PetscInt *cDof;
178: PetscSectionGetConstraintIndices(s, point, &cDof);
179: if (doInsert) {
180: for (i = 0; i < dim; ++i) {
181: if ((cInd < cDim) && (i == cDof[cInd])) {
182: if (doBC) array[i] = values[i]; /* Constrained update */
183: ++cInd;
184: continue;
185: }
186: if (doInterior) array[i] = values[i]; /* Unconstrained update */
187: }
188: } else {
189: for (i = 0; i < dim; ++i) {
190: if ((cInd < cDim) && (i == cDof[cInd])) {
191: if (doBC) array[i] += values[i]; /* Constrained update */
192: ++cInd;
193: continue;
194: }
195: if (doInterior) array[i] += values[i]; /* Unconstrained update */
196: }
197: }
198: } else {
199: /* TODO This is broken for add and constrained update */
200: const PetscInt *cDof;
201: PetscInt offset = 0;
202: PetscInt cOffset = 0;
203: PetscInt j = 0, field;
205: PetscSectionGetConstraintIndices(s, point, &cDof);
206: for (field = 0; field < s->numFields; ++field) {
207: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
208: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
209: const PetscInt sDim = dim - tDim;
210: PetscInt cInd = 0, i ,k;
212: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
213: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
214: if (doInterior) array[j] = values[k]; /* Unconstrained update */
215: }
216: offset += dim;
217: cOffset += dim - tDim;
218: }
219: }
220: }
221: VecRestoreArray(v, &baseArray);
222: return(0);
223: }