Actual source code: vsection.c
petsc-3.3-p7 2013-05-11
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc-private/vecimpl.h> /*I "petscvec.h" I*/
7: #if 0
8: /* Should I protect these for C++? */
11: PetscErrorCode PetscSectionGetDof(PetscUniformSection s, PetscInt point, PetscInt *numDof)
12: {
14: if ((point < s->pStart) || (point >= s->pEnd)) {
15: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
16: }
17: *numDof = s->numDof;
18: return(0);
19: }
23: PetscErrorCode PetscSectionGetOffset(PetscUniformSection s, PetscInt point, PetscInt *offset)
24: {
26: if ((point < s->pStart) || (point >= s->pEnd)) {
27: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
28: }
29: *offset = s->numDof*(point - s->pStart);
30: return(0);
31: }
32: #endif
36: /*@C
37: PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
39: Collective on MPI_Comm
41: Input Parameters:
42: + comm - the MPI communicator
43: - s - pointer to the section
45: Level: developer
47: Notes: Typical calling sequence
48: PetscSectionCreate(MPI_Comm,PetscSection *);
49: PetscSectionSetChart(PetscSection,low,high);
50: PetscSectionSetDof(PetscSection,point,numdof);
51: PetscSectionSetUp(PetscSection);
52: PetscSectionGetOffset(PetscSection,point,PetscInt *);
53: PetscSectionDestroy(PetscSection);
55: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
56: recommended they not be used in user codes unless you really gain something in their use.
58: Fortran Notes:
59: Not available from Fortran
61: .seealso: PetscSection, PetscSectionDestroy()
62: @*/
63: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
64: {
68: PetscNew(struct _n_PetscSection, s);
69: (*s)->atlasLayout.comm = comm;
70: (*s)->atlasLayout.pStart = -1;
71: (*s)->atlasLayout.pEnd = -1;
72: (*s)->atlasLayout.numDof = 1;
73: (*s)->atlasDof = PETSC_NULL;
74: (*s)->atlasOff = PETSC_NULL;
75: (*s)->bc = PETSC_NULL;
76: (*s)->bcIndices = PETSC_NULL;
77: (*s)->setup = PETSC_FALSE;
78: (*s)->numFields = 0;
79: (*s)->fieldNames = PETSC_NULL;
80: (*s)->field = PETSC_NULL;
81: return(0);
82: }
86: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
87: {
90: *numFields = s->numFields;
91: return(0);
92: }
96: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
97: {
98: PetscInt f;
102: if (numFields <= 0) {
103: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %d must be positive", numFields);
104: }
105: s->numFields = numFields;
106: PetscMalloc(s->numFields * sizeof(PetscInt), &s->numFieldComponents);
107: PetscMalloc(s->numFields * sizeof(char *), &s->fieldNames);
108: PetscMalloc(s->numFields * sizeof(PetscSection), &s->field);
109: for(f = 0; f < s->numFields; ++f) {
110: char name[64];
112: s->numFieldComponents[f] = 1;
113: PetscSectionCreate(s->atlasLayout.comm, &s->field[f]);
114: PetscSNPrintf(name, 64, "Field_%D", f);
115: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
116: }
117: return(0);
118: }
122: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
123: {
126: if ((field < 0) || (field >= s->numFields)) {
127: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
128: }
129: *fieldName = s->fieldNames[field];
130: return(0);
131: }
135: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
136: {
141: if ((field < 0) || (field >= s->numFields)) {
142: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
143: }
144: PetscFree(s->fieldNames[field]);
145: PetscStrallocpy(fieldName, (char **) &s->fieldNames[field]);
146: return(0);
147: }
151: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
152: {
155: if ((field < 0) || (field >= s->numFields)) {
156: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
157: }
158: *numComp = s->numFieldComponents[field];
159: return(0);
160: }
164: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
165: {
167: if ((field < 0) || (field >= s->numFields)) {
168: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
169: }
170: s->numFieldComponents[field] = numComp;
171: return(0);
172: }
176: PetscErrorCode PetscSectionCheckConstraints(PetscSection s)
177: {
181: if (!s->bc) {
182: PetscSectionCreate(s->atlasLayout.comm, &s->bc);
183: PetscSectionSetChart(s->bc, s->atlasLayout.pStart, s->atlasLayout.pEnd);
184: }
185: return(0);
186: }
190: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
191: {
193: if (pStart) {*pStart = s->atlasLayout.pStart;}
194: if (pEnd) {*pEnd = s->atlasLayout.pEnd;}
195: return(0);
196: }
200: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
201: {
202: PetscInt f;
206: s->atlasLayout.pStart = pStart;
207: s->atlasLayout.pEnd = pEnd;
208: PetscFree2(s->atlasDof, s->atlasOff);
209: PetscMalloc2((pEnd - pStart), PetscInt, &s->atlasDof, (pEnd - pStart), PetscInt, &s->atlasOff);
210: PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
211: for(f = 0; f < s->numFields; ++f) {
212: PetscSectionSetChart(s->field[f], pStart, pEnd);
213: }
214: return(0);
215: }
219: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
220: {
222: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) {
223: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
224: }
225: *numDof = s->atlasDof[point - s->atlasLayout.pStart];
226: return(0);
227: }
231: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
232: {
234: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) {
235: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
236: }
237: s->atlasDof[point - s->atlasLayout.pStart] = numDof;
238: return(0);
239: }
243: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
244: {
246: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) {
247: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
248: }
249: s->atlasDof[point - s->atlasLayout.pStart] += numDof;
250: return(0);
251: }
255: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
256: {
260: if ((field < 0) || (field >= s->numFields)) {
261: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
262: }
263: PetscSectionGetDof(s->field[field], point, numDof);
264: return(0);
265: }
269: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
270: {
274: if ((field < 0) || (field >= s->numFields)) {
275: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
276: }
277: PetscSectionSetDof(s->field[field], point, numDof);
278: return(0);
279: }
283: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
284: {
288: if (s->bc) {
289: PetscSectionGetDof(s->bc, point, numDof);
290: } else {
291: *numDof = 0;
292: }
293: return(0);
294: }
298: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
299: {
303: if (numDof) {
304: PetscSectionCheckConstraints(s);
305: PetscSectionSetDof(s->bc, point, numDof);
306: }
307: return(0);
308: }
312: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
313: {
317: if (numDof) {
318: PetscSectionCheckConstraints(s);
319: PetscSectionAddDof(s->bc, point, numDof);
320: }
321: return(0);
322: }
326: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
327: {
331: if ((field < 0) || (field >= s->numFields)) {
332: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
333: }
334: PetscSectionGetConstraintDof(s->field[field], point, numDof);
335: return(0);
336: }
340: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
341: {
345: if ((field < 0) || (field >= s->numFields)) {
346: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
347: }
348: PetscSectionSetConstraintDof(s->field[field], point, numDof);
349: return(0);
350: }
354: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
355: {
359: if (s->bc) {
360: const PetscInt last = (s->bc->atlasLayout.pEnd-s->bc->atlasLayout.pStart) - 1;
362: PetscSectionSetUp(s->bc);
363: PetscMalloc((s->bc->atlasOff[last] + s->bc->atlasDof[last]) * sizeof(PetscInt), &s->bcIndices);
364: }
365: return(0);
366: }
370: PetscErrorCode PetscSectionSetUp(PetscSection s)
371: {
372: PetscInt offset = 0, p, f;
376: if (s->setup) {return(0);}
377: s->setup = PETSC_TRUE;
378: for(p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
379: s->atlasOff[p] = offset;
380: offset += s->atlasDof[p];
381: }
382: PetscSectionSetUpBC(s);
383: /* Assume that all fields have the same chart */
384: for(p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
385: PetscInt off = s->atlasOff[p];
387: for(f = 0; f < s->numFields; ++f) {
388: PetscSection sf = s->field[f];
390: sf->atlasOff[p] = off;
391: off += sf->atlasDof[p];
392: }
393: }
394: for(f = 0; f < s->numFields; ++f) {
395: PetscSectionSetUpBC(s->field[f]);
396: }
397: return(0);
398: }
402: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
403: {
404: PetscInt p, n = 0;
407: for(p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
408: n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
409: }
410: *size = n;
411: return(0);
412: }
416: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
417: {
418: PetscInt p, n = 0;
421: for(p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
422: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
423: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
424: }
425: *size = n;
426: return(0);
427: }
431: /*
432: This gives negative offsets to points not owned by this process
433: */
434: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscSection *gsection)
435: {
436: PetscInt *neg;
437: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots;
441: PetscSectionCreate(s->atlasLayout.comm, gsection);
442: PetscSectionGetChart(s, &pStart, &pEnd);
443: PetscSectionSetChart(*gsection, pStart, pEnd);
444: PetscMalloc((pEnd - pStart) * sizeof(PetscInt), &neg);
445: /* Mark ghost points with negative dof */
446: for(p = pStart; p < pEnd; ++p) {
447: PetscSectionGetDof(s, p, &dof);
448: PetscSectionSetDof(*gsection, p, dof);
449: PetscSectionGetConstraintDof(s, p, &cdof);
450: if (cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
451: neg[p-pStart] = -(dof+1);
452: }
453: PetscSectionSetUpBC(*gsection);
454: PetscSFGetGraph(sf, &nroots, PETSC_NULL, PETSC_NULL, PETSC_NULL);
455: if (nroots >= 0) {
456: PetscSFBcastBegin(sf, MPIU_INT, &neg[-pStart], &(*gsection)->atlasDof[-pStart]);
457: PetscSFBcastEnd(sf, MPIU_INT, &neg[-pStart], &(*gsection)->atlasDof[-pStart]);
458: }
459: /* Calculate new sizes, get proccess offset, and calculate point offsets */
460: for(p = 0, off = 0; p < pEnd-pStart; ++p) {
461: cdof = s->bc ? s->bc->atlasDof[p] : 0;
462: (*gsection)->atlasOff[p] = off;
463: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
464: }
465: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, s->atlasLayout.comm);
466: globalOff -= off;
467: for(p = 0, off = 0; p < pEnd-pStart; ++p) {
468: (*gsection)->atlasOff[p] += globalOff;
469: neg[p] = -((*gsection)->atlasOff[p]+1);
470: }
471: /* Put in negative offsets for ghost points */
472: if (nroots >= 0) {
473: PetscSFBcastBegin(sf, MPIU_INT, &neg[-pStart], &(*gsection)->atlasOff[-pStart]);
474: PetscSFBcastEnd(sf, MPIU_INT, &neg[-pStart], &(*gsection)->atlasOff[-pStart]);
475: }
476: PetscFree(neg);
477: return(0);
478: }
482: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
483: {
484: PetscInt pStart, pEnd, p, localSize = 0;
488: PetscSectionGetChart(s, &pStart, &pEnd);
489: for(p = pStart; p < pEnd; ++p) {
490: PetscInt dof;
492: PetscSectionGetDof(s, p, &dof);
493: if (dof > 0) {++localSize;}
494: }
495: PetscLayoutCreate(comm, layout);
496: PetscLayoutSetLocalSize(*layout, localSize);
497: PetscLayoutSetBlockSize(*layout, 1);
498: PetscLayoutSetUp(*layout);
499: return(0);
500: }
504: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
505: {
506: PetscInt pStart, pEnd, p, localSize = 0;
510: PetscSectionGetChart(s, &pStart, &pEnd);
511: for(p = pStart; p < pEnd; ++p) {
512: PetscInt dof;
514: PetscSectionGetDof(s, p, &dof);
515: if (dof > 0) {localSize += dof;}
516: }
517: PetscLayoutCreate(comm, layout);
518: PetscLayoutSetLocalSize(*layout, localSize);
519: PetscLayoutSetBlockSize(*layout, 1);
520: PetscLayoutSetUp(*layout);
521: return(0);
522: }
526: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
527: {
529: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) {
530: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
531: }
532: *offset = s->atlasOff[point - s->atlasLayout.pStart];
533: return(0);
534: }
538: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
539: {
543: if ((field < 0) || (field >= s->numFields)) {
544: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
545: }
546: PetscSectionGetOffset(s->field[field], point, offset);
547: return(0);
548: }
552: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
553: {
554: PetscInt os = s->atlasOff[0], oe = s->atlasOff[0], pStart, pEnd, p;
558: PetscSectionGetChart(s, &pStart, &pEnd);
559: for(p = pStart; p < pEnd; ++p) {
560: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
562: if (off >= 0) {
563: os = PetscMin(os, off);
564: oe = PetscMax(oe, off+dof);
565: }
566: }
567: if (start) *start = os;
568: if (end) *end = oe;
569: return(0);
570: }
574: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
575: {
576: PetscInt p;
577: PetscMPIInt rank;
581: if (s->atlasLayout.numDof != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle %d dof in a uniform section", s->atlasLayout.numDof);
582: MPI_Comm_rank(((PetscObject) viewer)->comm, &rank);
583: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
584: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
585: for(p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
586: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
587: PetscInt b;
589: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d constrained", p+s->atlasLayout.pStart, s->atlasDof[p], s->atlasOff[p]);
590: for(b = 0; b < s->bc->atlasDof[p]; ++b) {
591: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
592: }
593: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
594: } else {
595: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d\n", p+s->atlasLayout.pStart, s->atlasDof[p], s->atlasOff[p]);
596: }
597: }
598: PetscViewerFlush(viewer);
599: return(0);
600: }
604: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
605: {
606: PetscBool isascii;
607: PetscInt f;
611: if (!viewer) {PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer);}
613: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
614: if (isascii) {
615: if (s->numFields) {
616: PetscViewerASCIIPrintf(viewer, "PetscSection with %d fields\n", s->numFields);
617: for(f = 0; f < s->numFields; ++f) {
618: PetscViewerASCIIPrintf(viewer, " field %d with %d components\n", f, s->numFieldComponents[f]);
619: PetscSectionView_ASCII(s->field[f], viewer);
620: }
621: } else {
622: PetscViewerASCIIPrintf(viewer, "PetscSection\n");
623: PetscSectionView_ASCII(s, viewer);
624: }
625: } else {
626: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported by this section object", ((PetscObject) viewer)->type_name);
627: }
628: return(0);
629: }
633: PetscErrorCode PetscSectionVecView_ASCII(PetscSection s, Vec v, PetscViewer viewer)
634: {
635: PetscScalar *array;
636: PetscInt p, i;
637: PetscMPIInt rank;
641: if (s->atlasLayout.numDof != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle %d dof in a uniform section", s->atlasLayout.numDof);
642: MPI_Comm_rank(((PetscObject) viewer)->comm, &rank);
643: VecGetArray(v, &array);
644: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
645: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
646: for(p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
647: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
648: PetscInt b;
650: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d", p+s->atlasLayout.pStart, s->atlasDof[p], s->atlasOff[p]);
651: for(i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
652: PetscScalar v = array[i];
653: #if defined(PETSC_USE_COMPLEX)
654: if (PetscImaginaryPart(v) > 0.0) {
655: PetscViewerASCIISynchronizedPrintf(viewer," %G + %G i", PetscRealPart(v), PetscImaginaryPart(v));
656: } else if (PetscImaginaryPart(v) < 0.0) {
657: PetscViewerASCIISynchronizedPrintf(viewer," %G - %G i", PetscRealPart(v),-PetscImaginaryPart(v));
658: } else {
659: PetscViewerASCIISynchronizedPrintf(viewer, " %G", PetscRealPart(v));
660: }
661: #else
662: PetscViewerASCIISynchronizedPrintf(viewer, " %G", v);
663: #endif
664: }
665: PetscViewerASCIISynchronizedPrintf(viewer, " constrained");
666: for(b = 0; b < s->bc->atlasDof[p]; ++b) {
667: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
668: }
669: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
670: } else {
671: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d", p+s->atlasLayout.pStart, s->atlasDof[p], s->atlasOff[p]);
672: for(i = s->atlasOff[p]; i < s->atlasOff[p]+s->atlasDof[p]; ++i) {
673: PetscScalar v = array[i];
674: #if defined(PETSC_USE_COMPLEX)
675: if (PetscImaginaryPart(v) > 0.0) {
676: PetscViewerASCIISynchronizedPrintf(viewer," %G + %G i", PetscRealPart(v), PetscImaginaryPart(v));
677: } else if (PetscImaginaryPart(v) < 0.0) {
678: PetscViewerASCIISynchronizedPrintf(viewer," %G - %G i", PetscRealPart(v),-PetscImaginaryPart(v));
679: } else {
680: PetscViewerASCIISynchronizedPrintf(viewer, " %G", PetscRealPart(v));
681: }
682: #else
683: PetscViewerASCIISynchronizedPrintf(viewer, " %G", v);
684: #endif
685: }
686: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
687: }
688: }
689: PetscViewerFlush(viewer);
690: VecRestoreArray(v, &array);
691: return(0);
692: }
696: PetscErrorCode PetscSectionVecView(PetscSection s, Vec v, PetscViewer viewer)
697: {
698: PetscBool isascii;
699: PetscInt f;
703: if (!viewer) {PetscViewerASCIIGetStdout(((PetscObject) v)->comm, &viewer);}
706: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
707: if (isascii) {
708: const char *name;
710: PetscObjectGetName((PetscObject) v, &name);
711: if (s->numFields) {
712: PetscViewerASCIIPrintf(viewer, "%s with %d fields\n", name, s->numFields);
713: for(f = 0; f < s->numFields; ++f) {
714: PetscViewerASCIIPrintf(viewer, " field %d with %d components\n", f, s->numFieldComponents[f]);
715: PetscSectionVecView_ASCII(s->field[f], v, viewer);
716: }
717: } else {
718: PetscViewerASCIIPrintf(viewer, "%s\n", name);
719: PetscSectionVecView_ASCII(s, v, viewer);
720: }
721: } else {
722: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported by this section object", ((PetscObject) viewer)->type_name);
723: }
724: return(0);
725: }
727: /*@C
728: PetscSectionDestroy - Frees a section object and frees its range if that exists.
730: Collective on MPI_Comm
732: Input Parameters:
733: . s - the PetscSection
735: Level: developer
737: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
738: recommended they not be used in user codes unless you really gain something in their use.
740: Fortran Notes:
741: Not available from Fortran
743: .seealso: PetscSection, PetscSectionCreate()
744: @*/
747: PetscErrorCode PetscSectionDestroy(PetscSection *s)
748: {
752: if (!*s) return(0);
753: if (!(*s)->refcnt--) {
754: PetscInt f;
756: PetscFree((*s)->numFieldComponents);
757: for(f = 0; f < (*s)->numFields; ++f) {
758: PetscSectionDestroy(&(*s)->field[f]);
759: PetscFree((*s)->fieldNames[f]);
760: }
761: PetscFree((*s)->fieldNames);
762: PetscFree((*s)->field);
763: PetscSectionDestroy(&(*s)->bc);
764: PetscFree((*s)->bcIndices);
765: PetscFree2((*s)->atlasDof, (*s)->atlasOff);
766: PetscFree((*s));
767: }
768: *s = PETSC_NULL;
769: return(0);
770: }
774: PetscErrorCode VecGetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar **values)
775: {
776: PetscScalar *baseArray;
777: const PetscInt p = point - s->atlasLayout.pStart;
781: VecGetArray(v, &baseArray);
782: *values = &baseArray[s->atlasOff[p]];
783: VecRestoreArray(v, &baseArray);
784: return(0);
785: }
789: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, PetscInt **values)
790: {
791: const PetscInt p = point - s->atlasLayout.pStart;
794: *values = &baseArray[s->atlasOff[p]];
795: return(0);
796: }
800: PetscErrorCode VecSetValuesSection(Vec v, PetscSection s, PetscInt point, PetscScalar values[], InsertMode mode)
801: {
802: PetscScalar *baseArray, *array;
803: const PetscBool doInsert = mode == INSERT_VALUES || mode == INSERT_ALL_VALUES ? PETSC_TRUE : PETSC_FALSE;
804: const PetscBool doBC = mode == INSERT_ALL_VALUES || mode == ADD_ALL_VALUES ? PETSC_TRUE : PETSC_FALSE;
805: const PetscInt p = point - s->atlasLayout.pStart;
806: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
807: PetscInt cDim = 0;
808: PetscErrorCode ierr;
811: PetscSectionGetConstraintDof(s, p, &cDim);
812: VecGetArray(v, &baseArray);
813: array = &baseArray[s->atlasOff[p]];
814: if (!cDim) {
815: if (orientation >= 0) {
816: const PetscInt dim = s->atlasDof[p];
817: PetscInt i;
819: if (doInsert) {
820: for(i = 0; i < dim; ++i) {
821: array[i] = values[i];
822: }
823: } else {
824: for(i = 0; i < dim; ++i) {
825: array[i] += values[i];
826: }
827: }
828: } else {
829: PetscInt offset = 0;
830: PetscInt j = -1, field, i;
832: for(field = 0; field < s->numFields; ++field) {
833: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
835: for(i = dim-1; i >= 0; --i) {
836: array[++j] = values[i+offset];
837: }
838: offset += dim;
839: }
840: }
841: } else {
842: if (orientation >= 0) {
843: const PetscInt dim = s->atlasDof[p];
844: PetscInt cInd = 0, i;
845: PetscInt *cDof;
847: PetscSectionGetConstraintIndices(s, point, &cDof);
848: if (doInsert) {
849: for(i = 0; i < dim; ++i) {
850: if ((cInd < cDim) && (i == cDof[cInd])) {
851: if (doBC) {array[i] = values[i];} /* Constrained update */
852: ++cInd;
853: continue;
854: }
855: array[i] = values[i]; /* Unconstrained update */
856: }
857: } else {
858: for(i = 0; i < dim; ++i) {
859: if ((cInd < cDim) && (i == cDof[cInd])) {
860: if (doBC) {array[i] += values[i];} /* Constrained update */
861: ++cInd;
862: continue;
863: }
864: array[i] += values[i]; /* Unconstrained update */
865: }
866: }
867: } else {
868: PetscInt *cDof;
869: PetscInt offset = 0;
870: PetscInt cOffset = 0;
871: PetscInt j = 0, field;
873: PetscSectionGetConstraintIndices(s, point, &cDof);
874: for(field = 0; field < s->numFields; ++field) {
875: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
876: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
877: const PetscInt sDim = dim - tDim;
878: PetscInt cInd = 0, i ,k;
880: for(i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
881: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
882: array[j] = values[k];
883: }
884: offset += dim;
885: cOffset += dim - tDim;
886: }
887: }
888: }
889: VecRestoreArray(v, &baseArray);
890: return(0);
891: }
896: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, PetscInt values[], InsertMode mode)
897: {
898: PetscInt *array;
899: const PetscInt p = point - s->atlasLayout.pStart;
900: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
901: PetscInt cDim = 0;
905: PetscSectionGetConstraintDof(s, p, &cDim);
906: array = &baseArray[s->atlasOff[p]];
907: if (!cDim) {
908: if (orientation >= 0) {
909: const PetscInt dim = s->atlasDof[p];
910: PetscInt i;
912: if (mode == INSERT_VALUES) {
913: for(i = 0; i < dim; ++i) {
914: array[i] = values[i];
915: }
916: } else {
917: for(i = 0; i < dim; ++i) {
918: array[i] += values[i];
919: }
920: }
921: } else {
922: PetscInt offset = 0;
923: PetscInt j = -1, field, i;
925: for(field = 0; field < s->numFields; ++field) {
926: const PetscInt dim = s->field[field]->atlasDof[p];
928: for(i = dim-1; i >= 0; --i) {
929: array[++j] = values[i+offset];
930: }
931: offset += dim;
932: }
933: }
934: } else {
935: if (orientation >= 0) {
936: const PetscInt dim = s->atlasDof[p];
937: PetscInt cInd = 0, i;
938: PetscInt *cDof;
940: PetscSectionGetConstraintIndices(s, point, &cDof);
941: if (mode == INSERT_VALUES) {
942: for(i = 0; i < dim; ++i) {
943: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
944: array[i] = values[i];
945: }
946: } else {
947: for(i = 0; i < dim; ++i) {
948: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
949: array[i] += values[i];
950: }
951: }
952: } else {
953: PetscInt *cDof;
954: PetscInt offset = 0;
955: PetscInt cOffset = 0;
956: PetscInt j = 0, field;
958: PetscSectionGetConstraintIndices(s, point, &cDof);
959: for(field = 0; field < s->numFields; ++field) {
960: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
961: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
962: const PetscInt sDim = dim - tDim;
963: PetscInt cInd = 0, i ,k;
965: for(i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
966: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
967: array[j] = values[k];
968: }
969: offset += dim;
970: cOffset += dim - tDim;
971: }
972: }
973: }
974: return(0);
975: }
979: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, PetscInt **indices)
980: {
984: if (s->bc) {
985: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
986: } else {
987: *indices = PETSC_NULL;
988: }
989: return(0);
990: }
994: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, PetscInt indices[])
995: {
999: if (s->bc) {
1000: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1001: }
1002: return(0);
1003: }
1007: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, PetscInt **indices)
1008: {
1012: if ((field < 0) || (field >= s->numFields)) {
1013: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1014: }
1015: PetscSectionGetConstraintIndices(s->field[field], point, indices);
1016: return(0);
1017: }
1021: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, PetscInt indices[])
1022: {
1026: if ((field < 0) || (field >= s->numFields)) {
1027: SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1028: }
1029: PetscSectionSetConstraintIndices(s->field[field], point, indices);
1030: return(0);
1031: }
1035: PetscErrorCode PetscSFConvertPartition(PetscSF sfPart, PetscSection partSection, IS partition, ISLocalToGlobalMapping *renumbering, PetscSF *sf)
1036: {
1037: MPI_Comm comm = ((PetscObject)sfPart)->comm;
1038: PetscSF sfPoints;
1039: PetscInt *partSizes,*partOffsets,p,i,numParts,numMyPoints,numPoints,count;
1040: const PetscInt *partArray;
1041: PetscSFNode *sendPoints;
1042: PetscMPIInt rank;
1046: MPI_Comm_rank(comm, &rank);
1048: /* Get the number of parts and sizes that I have to distribute */
1049: PetscSFGetGraph(sfPart,PETSC_NULL,&numParts,PETSC_NULL,PETSC_NULL);
1050: PetscMalloc2(numParts,PetscInt,&partSizes,numParts,PetscInt,&partOffsets);
1051: for (p=0,numPoints=0; p<numParts; p++) {
1052: PetscSectionGetDof(partSection, p, &partSizes[p]);
1053: numPoints += partSizes[p];
1054: }
1055: numMyPoints = 0;
1056: PetscSFFetchAndOpBegin(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1057: PetscSFFetchAndOpEnd(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1058: /* I will receive numMyPoints. I will send a total of numPoints, to be placed on remote procs at partOffsets */
1060: /* Create SF mapping locations (addressed through partition, as indexed leaves) to new owners (roots) */
1061: PetscMalloc(numPoints*sizeof(PetscSFNode),&sendPoints);
1062: for (p=0,count=0; p<numParts; p++) {
1063: for (i=0; i<partSizes[p]; i++) {
1064: sendPoints[count].rank = p;
1065: sendPoints[count].index = partOffsets[p]+i;
1066: count++;
1067: }
1068: }
1069: if (count != numPoints) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count %D should equal numPoints=%D",count,numPoints);
1070: PetscFree2(partSizes,partOffsets);
1071: ISGetIndices(partition,&partArray);
1072: PetscSFCreate(comm,&sfPoints);
1073: PetscSFSetGraph(sfPoints,numMyPoints,numPoints,partArray,PETSC_USE_POINTER,sendPoints,PETSC_OWN_POINTER);
1075: /* Invert SF so that the new owners are leaves and the locations indexed through partition are the roots */
1076: PetscSFCreateInverseSF(sfPoints,sf);
1077: PetscSFDestroy(&sfPoints);
1078: ISRestoreIndices(partition,&partArray);
1080: /* Create the new local-to-global mapping */
1081: ISLocalToGlobalMappingCreateSF(*sf,0,renumbering);
1082: return(0);
1083: }
1087: /*@C
1088: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
1090: Collective
1092: Input Parameters:
1093: + sf - The SF
1094: - rootSection - Section defined on root space
1096: Output Parameters:
1097: + remoteOffsets - root offsets in leaf storage, or PETSC_NULL
1098: - leafSection - Section defined on the leaf space
1100: Level: intermediate
1102: .seealso: PetscSFCreate()
1103: @*/
1104: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1105: {
1106: PetscSF embedSF;
1107: const PetscInt *ilocal, *indices;
1108: IS selected;
1109: PetscInt nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i;
1110: PetscErrorCode ierr;
1113: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1114: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1115: ISGetIndices(selected, &indices);
1116: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1117: ISRestoreIndices(selected, &indices);
1118: ISDestroy(&selected);
1119: PetscSFGetGraph(embedSF, PETSC_NULL, &nleaves, &ilocal, PETSC_NULL);
1120: if (ilocal) {
1121: for(i = 0; i < nleaves; ++i) {
1122: lpStart = PetscMin(lpStart, ilocal[i]);
1123: lpEnd = PetscMax(lpEnd, ilocal[i]);
1124: }
1125: } else {
1126: lpStart = 0;
1127: lpEnd = nleaves;
1128: }
1129: ++lpEnd;
1130: PetscSectionSetChart(leafSection, lpStart, lpEnd);
1131: /* Could fuse these at the cost of a copy and extra allocation */
1132: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1133: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1134: if (remoteOffsets) {
1135: PetscMalloc((lpEnd - lpStart) * sizeof(PetscInt), remoteOffsets);
1136: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1137: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1138: }
1139: PetscSFDestroy(&embedSF);
1140: PetscSectionSetUp(leafSection);
1141: return(0);
1142: }
1146: /*@C
1147: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
1149: Input Parameters:
1150: + sf - The SF
1151: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section), or PETSC_NULL
1152: - remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or PETSC_NULL
1154: Output Parameters:
1155: + leafSection - Data layout of local points for incoming data (this is the distributed section)
1156: - sectionSF - The new SF
1158: Note: Either rootSection or remoteOffsets can be specified
1160: Level: intermediate
1162: .seealso: PetscSFCreate()
1163: @*/
1164: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1165: {
1166: MPI_Comm comm = ((PetscObject) sf)->comm;
1167: const PetscInt *localPoints;
1168: const PetscSFNode *remotePoints;
1169: PetscInt lpStart, lpEnd;
1170: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
1171: PetscInt *localIndices;
1172: PetscSFNode *remoteIndices;
1173: PetscInt i, ind;
1174: PetscErrorCode ierr;
1177: PetscSFCreate(comm, sectionSF);
1178: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1179: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
1180: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
1181: if (numRoots < 0) {return(0);}
1182: for(i = 0; i < numPoints; ++i) {
1183: PetscInt localPoint = localPoints ? localPoints[i] : i;
1184: PetscInt dof;
1186: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
1187: PetscSectionGetDof(leafSection, localPoint, &dof);
1188: numIndices += dof;
1189: }
1190: }
1191: PetscMalloc(numIndices * sizeof(PetscInt), &localIndices);
1192: PetscMalloc(numIndices * sizeof(PetscSFNode), &remoteIndices);
1193: /* Get offsets for remote data */
1194: if (!remoteOffsets) {
1195: PetscSF embedSF;
1196: const PetscInt *indices;
1197: IS selected;
1198: PetscInt rpStart, rpEnd, isSize;
1200: PetscMalloc((lpEnd - lpStart) * sizeof(PetscInt), &remoteOffsets);
1201: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1202: isSize = PetscMin(numRoots, rpEnd - rpStart);
1203: ISCreateStride(PETSC_COMM_SELF, isSize, rpStart, 1, &selected);
1204: ISGetIndices(selected, &indices);
1205: #if 0
1206: PetscSFCreateEmbeddedSF(sf, isSize, indices, &embedSF);
1207: #else
1208: embedSF = sf;
1209: PetscObjectReference((PetscObject) embedSF);
1210: #endif
1211: ISRestoreIndices(selected, &indices);
1212: ISDestroy(&selected);
1213: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &remoteOffsets[-lpStart]);
1214: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &remoteOffsets[-lpStart]);
1215: PetscSFDestroy(&embedSF);
1216: }
1217: /* Create new index graph */
1218: for(i = 0, ind = 0; i < numPoints; ++i) {
1219: PetscInt localPoint = localPoints ? localPoints[i] : i;
1220: PetscInt rank = remotePoints[i].rank;
1222: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
1223: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
1224: PetscInt loff, dof, d;
1226: PetscSectionGetOffset(leafSection, localPoint, &loff);
1227: PetscSectionGetDof(leafSection, localPoint, &dof);
1228: for(d = 0; d < dof; ++d, ++ind) {
1229: localIndices[ind] = loff+d;
1230: remoteIndices[ind].rank = rank;
1231: remoteIndices[ind].index = remoteOffset+d;
1232: }
1233: }
1234: }
1235: PetscFree(remoteOffsets);
1236: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
1237: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
1238: return(0);
1239: }