Actual source code: vsectionis.c
petsc-3.4.5 2014-06-29
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc-private/isimpl.h> /*I "petscvec.h" I*/
6: #include <petscsf.h>
7: #include <petscviewer.h>
9: PetscClassId PETSC_SECTION_CLASSID;
13: /*@
14: PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
16: Collective on MPI_Comm
18: Input Parameters:
19: + comm - the MPI communicator
20: - s - pointer to the section
22: Level: developer
24: Notes: Typical calling sequence
25: PetscSectionCreate(MPI_Comm,PetscSection *);
26: PetscSectionSetChart(PetscSection,low,high);
27: PetscSectionSetDof(PetscSection,point,numdof);
28: PetscSectionSetUp(PetscSection);
29: PetscSectionGetOffset(PetscSection,point,PetscInt *);
30: PetscSectionDestroy(PetscSection);
32: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
33: recommended they not be used in user codes unless you really gain something in their use.
35: .seealso: PetscSection, PetscSectionDestroy()
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
43: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
44: ISInitializePackage();
45: #endif
47: PetscHeaderCreate(*s,_p_PetscSection,int,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
49: (*s)->atlasLayout.comm = comm;
50: (*s)->atlasLayout.pStart = -1;
51: (*s)->atlasLayout.pEnd = -1;
52: (*s)->atlasLayout.numDof = 1;
53: (*s)->maxDof = 0;
54: (*s)->atlasDof = NULL;
55: (*s)->atlasOff = NULL;
56: (*s)->bc = NULL;
57: (*s)->bcIndices = NULL;
58: (*s)->setup = PETSC_FALSE;
59: (*s)->numFields = 0;
60: (*s)->fieldNames = NULL;
61: (*s)->field = NULL;
62: return(0);
63: }
67: /*@
68: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
70: Collective on MPI_Comm
72: Input Parameter:
73: . section - the PetscSection
75: Output Parameter:
76: . newSection - the copy
78: Level: developer
80: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
81: @*/
82: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
83: {
84: PetscInt numFields, f, pStart, pEnd, p;
88: PetscSectionCreate(section->atlasLayout.comm, newSection);
89: PetscSectionGetNumFields(section, &numFields);
90: if (numFields) {PetscSectionSetNumFields(*newSection, numFields);}
91: for (f = 0; f < numFields; ++f) {
92: const char *name = NULL;
93: PetscInt numComp = 0;
95: PetscSectionGetFieldName(section, f, &name);
96: PetscSectionSetFieldName(*newSection, f, name);
97: PetscSectionGetFieldComponents(section, f, &numComp);
98: PetscSectionSetFieldComponents(*newSection, f, numComp);
99: }
100: PetscSectionGetChart(section, &pStart, &pEnd);
101: PetscSectionSetChart(*newSection, pStart, pEnd);
102: for (p = pStart; p < pEnd; ++p) {
103: PetscInt dof, cdof, fcdof = 0;
105: PetscSectionGetDof(section, p, &dof);
106: PetscSectionSetDof(*newSection, p, dof);
107: PetscSectionGetConstraintDof(section, p, &cdof);
108: if (cdof) {PetscSectionSetConstraintDof(*newSection, p, cdof);}
109: for (f = 0; f < numFields; ++f) {
110: PetscSectionGetFieldDof(section, p, f, &dof);
111: PetscSectionSetFieldDof(*newSection, p, f, dof);
112: if (cdof) {
113: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
114: if (fcdof) {PetscSectionSetFieldConstraintDof(*newSection, p, f, fcdof);}
115: }
116: }
117: }
118: PetscSectionSetUp(*newSection);
119: for (p = pStart; p < pEnd; ++p) {
120: PetscInt off, cdof, fcdof = 0;
121: const PetscInt *cInd;
123: /* Must set offsets in case they do not agree with the prefix sums */
124: PetscSectionGetOffset(section, p, &off);
125: PetscSectionSetOffset(*newSection, p, off);
126: PetscSectionGetConstraintDof(section, p, &cdof);
127: if (cdof) {
128: PetscSectionGetConstraintIndices(section, p, &cInd);
129: PetscSectionSetConstraintIndices(*newSection, p, cInd);
130: for (f = 0; f < numFields; ++f) {
131: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
132: if (fcdof) {
133: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
134: PetscSectionSetFieldConstraintIndices(*newSection, p, f, cInd);
135: }
136: }
137: }
138: }
139: return(0);
140: }
144: /*@
145: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
147: Not collective
149: Input Parameter:
150: . s - the PetscSection
152: Output Parameter:
153: . numFields - the number of fields defined, or 0 if none were defined
155: Level: intermediate
157: .seealso: PetscSectionSetNumFields()
158: @*/
159: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
160: {
163: *numFields = s->numFields;
164: return(0);
165: }
169: /*@
170: PetscSectionSetNumFields - Sets the number of fields.
172: Not collective
174: Input Parameters:
175: + s - the PetscSection
176: - numFields - the number of fields
178: Level: intermediate
180: .seealso: PetscSectionGetNumFields()
181: @*/
182: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
183: {
184: PetscInt f;
188: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %d must be positive", numFields);
189: s->numFields = numFields;
191: PetscMalloc(s->numFields * sizeof(PetscInt), &s->numFieldComponents);
192: PetscMalloc(s->numFields * sizeof(char*), &s->fieldNames);
193: PetscMalloc(s->numFields * sizeof(PetscSection), &s->field);
194: for (f = 0; f < s->numFields; ++f) {
195: char name[64];
197: s->numFieldComponents[f] = 1;
199: PetscSectionCreate(s->atlasLayout.comm, &s->field[f]);
200: PetscSNPrintf(name, 64, "Field_%D", f);
201: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
202: }
203: return(0);
204: }
208: /*@C
209: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
211: Not Collective
213: Input Parameters:
214: + s - the PetscSection
215: - field - the field number
217: Output Parameter:
218: . fieldName - the field name
220: Level: developer
222: .seealso: PetscSectionSetFieldName()
223: @*/
224: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
225: {
228: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
229: *fieldName = s->fieldNames[field];
230: return(0);
231: }
235: /*@C
236: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
238: Not Collective
240: Input Parameters:
241: + s - the PetscSection
242: . field - the field number
243: - fieldName - the field name
245: Level: developer
247: .seealso: PetscSectionGetFieldName()
248: @*/
249: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
250: {
255: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
256: PetscFree(s->fieldNames[field]);
257: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
258: return(0);
259: }
263: /*@
264: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
266: Not collective
268: Input Parameters:
269: + s - the PetscSection
270: - field - the field number
272: Output Parameter:
273: . numComp - the number of field components
275: Level: intermediate
277: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
278: @*/
279: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
280: {
283: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
284: *numComp = s->numFieldComponents[field];
285: return(0);
286: }
290: /*@
291: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
293: Not collective
295: Input Parameters:
296: + s - the PetscSection
297: . field - the field number
298: - numComp - the number of field components
300: Level: intermediate
302: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
303: @*/
304: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
305: {
307: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
308: s->numFieldComponents[field] = numComp;
309: return(0);
310: }
314: PetscErrorCode PetscSectionCheckConstraints(PetscSection s)
315: {
319: if (!s->bc) {
320: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
321: PetscSectionSetChart(s->bc, s->atlasLayout.pStart, s->atlasLayout.pEnd);
322: }
323: return(0);
324: }
328: /*@
329: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points in the lie.
331: Not collective
333: Input Parameter:
334: . s - the PetscSection
336: Output Parameters:
337: + pStart - the first point
338: - pEnd - one past the last point
340: Level: intermediate
342: .seealso: PetscSectionSetChart(), PetscSectionCreate()
343: @*/
344: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
345: {
347: if (pStart) *pStart = s->atlasLayout.pStart;
348: if (pEnd) *pEnd = s->atlasLayout.pEnd;
349: return(0);
350: }
354: /*@
355: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points in the lie.
357: Not collective
359: Input Parameters:
360: + s - the PetscSection
361: . pStart - the first point
362: - pEnd - one past the last point
364: Level: intermediate
366: .seealso: PetscSectionSetChart(), PetscSectionCreate()
367: @*/
368: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
369: {
370: PetscInt f;
374: /* Cannot Reset() because it destroys field information */
375: s->setup = PETSC_FALSE;
376: PetscSectionDestroy(&s->bc);
377: PetscFree(s->bcIndices);
378: PetscFree2(s->atlasDof, s->atlasOff);
380: s->atlasLayout.pStart = pStart;
381: s->atlasLayout.pEnd = pEnd;
382: PetscMalloc2((pEnd - pStart), PetscInt, &s->atlasDof, (pEnd - pStart), PetscInt, &s->atlasOff);
383: PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
384: for (f = 0; f < s->numFields; ++f) {
385: PetscSectionSetChart(s->field[f], pStart, pEnd);
386: }
387: return(0);
388: }
392: /*@
393: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
395: Not collective
397: Input Parameters:
398: + s - the PetscSection
399: - point - the point
401: Output Parameter:
402: . numDof - the number of dof
404: Level: intermediate
406: .seealso: PetscSectionSetDof(), PetscSectionCreate()
407: @*/
408: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
409: {
411: #if defined(PETSC_USE_DEBUG)
412: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
413: #endif
414: *numDof = s->atlasDof[point - s->atlasLayout.pStart];
415: return(0);
416: }
420: /*@
421: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
423: Not collective
425: Input Parameters:
426: + s - the PetscSection
427: . point - the point
428: - numDof - the number of dof
430: Level: intermediate
432: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
433: @*/
434: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
435: {
437: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
438: s->atlasDof[point - s->atlasLayout.pStart] = numDof;
439: return(0);
440: }
444: /*@
445: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
447: Not collective
449: Input Parameters:
450: + s - the PetscSection
451: . point - the point
452: - numDof - the number of additional dof
454: Level: intermediate
456: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
457: @*/
458: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
459: {
461: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
462: s->atlasDof[point - s->atlasLayout.pStart] += numDof;
463: return(0);
464: }
468: /*@
469: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
471: Not collective
473: Input Parameters:
474: + s - the PetscSection
475: . point - the point
476: - field - the field
478: Output Parameter:
479: . numDof - the number of dof
481: Level: intermediate
483: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
484: @*/
485: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
486: {
490: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
491: PetscSectionGetDof(s->field[field], point, numDof);
492: return(0);
493: }
497: /*@
498: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
500: Not collective
502: Input Parameters:
503: + s - the PetscSection
504: . point - the point
505: . field - the field
506: - numDof - the number of dof
508: Level: intermediate
510: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
511: @*/
512: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
513: {
517: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
518: PetscSectionSetDof(s->field[field], point, numDof);
519: return(0);
520: }
524: /*@
525: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
527: Not collective
529: Input Parameters:
530: + s - the PetscSection
531: . point - the point
532: . field - the field
533: - numDof - the number of dof
535: Level: intermediate
537: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
538: @*/
539: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
540: {
544: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
545: PetscSectionAddDof(s->field[field], point, numDof);
546: return(0);
547: }
551: /*@
552: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
554: Not collective
556: Input Parameters:
557: + s - the PetscSection
558: - point - the point
560: Output Parameter:
561: . numDof - the number of dof which are fixed by constraints
563: Level: intermediate
565: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
566: @*/
567: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
568: {
572: if (s->bc) {
573: PetscSectionGetDof(s->bc, point, numDof);
574: } else *numDof = 0;
575: return(0);
576: }
580: /*@
581: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
583: Not collective
585: Input Parameters:
586: + s - the PetscSection
587: . point - the point
588: - numDof - the number of dof which are fixed by constraints
590: Level: intermediate
592: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
593: @*/
594: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
595: {
599: if (numDof) {
600: PetscSectionCheckConstraints(s);
601: PetscSectionSetDof(s->bc, point, numDof);
602: }
603: return(0);
604: }
608: /*@
609: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
611: Not collective
613: Input Parameters:
614: + s - the PetscSection
615: . point - the point
616: - numDof - the number of additional dof which are fixed by constraints
618: Level: intermediate
620: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
621: @*/
622: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
623: {
627: if (numDof) {
628: PetscSectionCheckConstraints(s);
629: PetscSectionAddDof(s->bc, point, numDof);
630: }
631: return(0);
632: }
636: /*@
637: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
639: Not collective
641: Input Parameters:
642: + s - the PetscSection
643: . point - the point
644: - field - the field
646: Output Parameter:
647: . numDof - the number of dof which are fixed by constraints
649: Level: intermediate
651: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
652: @*/
653: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
654: {
658: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
659: PetscSectionGetConstraintDof(s->field[field], point, numDof);
660: return(0);
661: }
665: /*@
666: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
668: Not collective
670: Input Parameters:
671: + s - the PetscSection
672: . point - the point
673: . field - the field
674: - numDof - the number of dof which are fixed by constraints
676: Level: intermediate
678: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
679: @*/
680: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
681: {
685: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
686: PetscSectionSetConstraintDof(s->field[field], point, numDof);
687: return(0);
688: }
692: /*@
693: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
695: Not collective
697: Input Parameters:
698: + s - the PetscSection
699: . point - the point
700: . field - the field
701: - numDof - the number of additional dof which are fixed by constraints
703: Level: intermediate
705: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
706: @*/
707: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
708: {
712: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
713: PetscSectionAddConstraintDof(s->field[field], point, numDof);
714: return(0);
715: }
719: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
720: {
724: if (s->bc) {
725: const PetscInt last = (s->bc->atlasLayout.pEnd-s->bc->atlasLayout.pStart) - 1;
727: PetscSectionSetUp(s->bc);
728: PetscMalloc((s->bc->atlasOff[last] + s->bc->atlasDof[last]) * sizeof(PetscInt), &s->bcIndices);
729: }
730: return(0);
731: }
735: /*@
736: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
738: Not collective
740: Input Parameter:
741: . s - the PetscSection
743: Level: intermediate
745: .seealso: PetscSectionCreate()
746: @*/
747: PetscErrorCode PetscSectionSetUp(PetscSection s)
748: {
749: PetscInt offset = 0, p, f;
753: if (s->setup) return(0);
754: s->setup = PETSC_TRUE;
755: for (p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
756: s->atlasOff[p] = offset;
757: offset += s->atlasDof[p];
758: s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
759: }
760: PetscSectionSetUpBC(s);
761: /* Assume that all fields have the same chart */
762: for (p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
763: PetscInt off = s->atlasOff[p];
765: for (f = 0; f < s->numFields; ++f) {
766: PetscSection sf = s->field[f];
768: sf->atlasOff[p] = off;
769: off += sf->atlasDof[p];
770: }
771: }
772: for (f = 0; f < s->numFields; ++f) {
773: PetscSectionSetUpBC(s->field[f]);
774: }
775: return(0);
776: }
780: /*@
781: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
783: Not collective
785: Input Parameters:
786: . s - the PetscSection
788: Output Parameter:
789: . maxDof - the maximum dof
791: Level: intermediate
793: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
794: @*/
795: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
796: {
798: *maxDof = s->maxDof;
799: return(0);
800: }
804: /*@
805: PetscSectionGetStorageSize - Return the size of an arary or local Vec capable of holding all the degrees of freedom.
807: Not collective
809: Input Parameters:
810: + s - the PetscSection
811: - point - the point
813: Output Parameter:
814: . size - the size of an array which can hold all the dofs
816: Level: intermediate
818: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
819: @*/
820: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
821: {
822: PetscInt p, n = 0;
825: for (p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
826: *size = n;
827: return(0);
828: }
832: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
833: {
834: PetscInt p, n = 0;
837: for (p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
838: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
839: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
840: }
841: *size = n;
842: return(0);
843: }
847: /*@
848: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
849: the local section and an SF describing the section point overlap.
851: Input Parameters:
852: + s - The PetscSection for the local field layout
853: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
854: - includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
856: Output Parameter:
857: . gsection - The PetscSection for the global field layout
859: Note: This gives negative sizes and offsets to points not owned by this process
861: Level: developer
863: .seealso: PetscSectionCreate()
864: @*/
865: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscSection *gsection)
866: {
867: PetscInt *neg = NULL, *recv = NULL;
868: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
872: PetscSectionCreate(s->atlasLayout.comm, gsection);
873: PetscSectionGetChart(s, &pStart, &pEnd);
874: PetscSectionSetChart(*gsection, pStart, pEnd);
875: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
876: nlocal = nroots; /* The local/leaf space matches global/root space */
877: /* Must allocate for all points visible to SF, which may be more than this section */
878: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
879: PetscMalloc2(nroots,PetscInt,&neg,nlocal,PetscInt,&recv);
880: PetscMemzero(neg,nroots*sizeof(PetscInt));
881: }
882: /* Mark all local points with negative dof */
883: for (p = pStart; p < pEnd; ++p) {
884: PetscSectionGetDof(s, p, &dof);
885: PetscSectionSetDof(*gsection, p, dof);
886: PetscSectionGetConstraintDof(s, p, &cdof);
887: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
888: if (neg) neg[p] = -(dof+1);
889: }
890: PetscSectionSetUpBC(*gsection);
891: if (nroots >= 0) {
892: PetscMemzero(recv,nlocal*sizeof(PetscInt));
893: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
894: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
895: for (p = pStart; p < pEnd; ++p) {
896: if (recv[p] < 0) (*gsection)->atlasDof[p-pStart] = recv[p];
897: }
898: }
899: /* Calculate new sizes, get proccess offset, and calculate point offsets */
900: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
901: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
902: (*gsection)->atlasOff[p] = off;
903: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
904: }
905: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, s->atlasLayout.comm);
906: globalOff -= off;
907: for (p = pStart, off = 0; p < pEnd; ++p) {
908: (*gsection)->atlasOff[p-pStart] += globalOff;
909: if (neg) neg[p] = -((*gsection)->atlasOff[p-pStart]+1);
910: }
911: /* Put in negative offsets for ghost points */
912: if (nroots >= 0) {
913: PetscMemzero(recv,nlocal*sizeof(PetscInt));
914: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
915: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
916: for (p = pStart; p < pEnd; ++p) {
917: if (recv[p] < 0) (*gsection)->atlasOff[p-pStart] = recv[p];
918: }
919: }
920: PetscFree2(neg,recv);
921: return(0);
922: }
926: /*@
927: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
928: the local section and an SF describing the section point overlap.
930: Input Parameters:
931: + s - The PetscSection for the local field layout
932: . sf - The SF describing parallel layout of the section points
933: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
934: . numExcludes - The number of exclusion ranges
935: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
937: Output Parameter:
938: . gsection - The PetscSection for the global field layout
940: Note: This gives negative sizes and offsets to points not owned by this process
942: Level: developer
944: .seealso: PetscSectionCreate()
945: @*/
946: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
947: {
948: PetscInt *neg;
949: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
953: PetscSectionCreate(s->atlasLayout.comm, gsection);
954: PetscSectionGetChart(s, &pStart, &pEnd);
955: PetscSectionSetChart(*gsection, pStart, pEnd);
956: PetscMalloc((pEnd - pStart) * sizeof(PetscInt), &neg);
957: /* Mark ghost points with negative dof */
958: for (p = pStart; p < pEnd; ++p) {
959: for (e = 0; e < numExcludes; ++e) {
960: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
961: PetscSectionSetDof(*gsection, p, 0);
962: break;
963: }
964: }
965: if (e < numExcludes) continue;
966: PetscSectionGetDof(s, p, &dof);
967: PetscSectionSetDof(*gsection, p, dof);
968: PetscSectionGetConstraintDof(s, p, &cdof);
969: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
970: neg[p-pStart] = -(dof+1);
971: }
972: PetscSectionSetUpBC(*gsection);
973: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
974: if (nroots >= 0) {
975: if (nroots > pEnd - pStart) {
976: PetscInt *tmpDof;
977: /* Help Jed: HAVE TO MAKE A BUFFER HERE THE SIZE OF THE COMPLETE SPACE AND THEN COPY INTO THE atlasDof FOR THIS SECTION */
978: PetscMalloc(nroots * sizeof(PetscInt), &tmpDof);
979: PetscSFBcastBegin(sf, MPIU_INT, &neg[-pStart], tmpDof);
980: PetscSFBcastEnd(sf, MPIU_INT, &neg[-pStart], tmpDof);
981: for (p = pStart; p < pEnd; ++p) {
982: if (tmpDof[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpDof[p];
983: }
984: PetscFree(tmpDof);
985: } else {
986: PetscSFBcastBegin(sf, MPIU_INT, &neg[-pStart], &(*gsection)->atlasDof[-pStart]);
987: PetscSFBcastEnd(sf, MPIU_INT, &neg[-pStart], &(*gsection)->atlasDof[-pStart]);
988: }
989: }
990: /* Calculate new sizes, get proccess offset, and calculate point offsets */
991: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
992: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[p] : 0;
993: (*gsection)->atlasOff[p] = off;
994: off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
995: }
996: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, s->atlasLayout.comm);
997: globalOff -= off;
998: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
999: (*gsection)->atlasOff[p] += globalOff;
1000: neg[p] = -((*gsection)->atlasOff[p]+1);
1001: }
1002: /* Put in negative offsets for ghost points */
1003: if (nroots >= 0) {
1004: if (nroots > pEnd - pStart) {
1005: PetscInt *tmpOff;
1006: /* Help Jed: HAVE TO MAKE A BUFFER HERE THE SIZE OF THE COMPLETE SPACE AND THEN COPY INTO THE atlasDof FOR THIS SECTION */
1007: PetscMalloc(nroots * sizeof(PetscInt), &tmpOff);
1008: PetscSFBcastBegin(sf, MPIU_INT, &neg[-pStart], tmpOff);
1009: PetscSFBcastEnd(sf, MPIU_INT, &neg[-pStart], tmpOff);
1010: for (p = pStart; p < pEnd; ++p) {
1011: if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];
1012: }
1013: PetscFree(tmpOff);
1014: } else {
1015: PetscSFBcastBegin(sf, MPIU_INT, &neg[-pStart], &(*gsection)->atlasOff[-pStart]);
1016: PetscSFBcastEnd(sf, MPIU_INT, &neg[-pStart], &(*gsection)->atlasOff[-pStart]);
1017: }
1018: }
1019: PetscFree(neg);
1020: return(0);
1021: }
1025: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1026: {
1027: PetscInt pStart, pEnd, p, localSize = 0;
1031: PetscSectionGetChart(s, &pStart, &pEnd);
1032: for (p = pStart; p < pEnd; ++p) {
1033: PetscInt dof;
1035: PetscSectionGetDof(s, p, &dof);
1036: if (dof > 0) ++localSize;
1037: }
1038: PetscLayoutCreate(comm, layout);
1039: PetscLayoutSetLocalSize(*layout, localSize);
1040: PetscLayoutSetBlockSize(*layout, 1);
1041: PetscLayoutSetUp(*layout);
1042: return(0);
1043: }
1047: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1048: {
1049: PetscInt pStart, pEnd, p, localSize = 0;
1053: PetscSectionGetChart(s, &pStart, &pEnd);
1054: for (p = pStart; p < pEnd; ++p) {
1055: PetscInt dof,cdof;
1057: PetscSectionGetDof(s, p, &dof);
1058: PetscSectionGetConstraintDof(s, p, &cdof);
1059: if (dof-cdof > 0) localSize += dof-cdof;
1060: }
1061: PetscLayoutCreate(comm, layout);
1062: PetscLayoutSetLocalSize(*layout, localSize);
1063: PetscLayoutSetBlockSize(*layout, 1);
1064: PetscLayoutSetUp(*layout);
1065: return(0);
1066: }
1070: /*@
1071: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1073: Not collective
1075: Input Parameters:
1076: + s - the PetscSection
1077: - point - the point
1079: Output Parameter:
1080: . offset - the offset
1082: Level: intermediate
1084: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1085: @*/
1086: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1087: {
1089: #if defined(PETSC_USE_DEBUG)
1090: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
1091: #endif
1092: *offset = s->atlasOff[point - s->atlasLayout.pStart];
1093: return(0);
1094: }
1098: /*@
1099: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1101: Not collective
1103: Input Parameters:
1104: + s - the PetscSection
1105: . point - the point
1106: - offset - the offset
1108: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1110: Level: intermediate
1112: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1113: @*/
1114: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1115: {
1117: if ((point < s->atlasLayout.pStart) || (point >= s->atlasLayout.pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->atlasLayout.pStart, s->atlasLayout.pEnd);
1118: s->atlasOff[point - s->atlasLayout.pStart] = offset;
1119: return(0);
1120: }
1124: /*@
1125: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1127: Not collective
1129: Input Parameters:
1130: + s - the PetscSection
1131: . point - the point
1132: - field - the field
1134: Output Parameter:
1135: . offset - the offset
1137: Level: intermediate
1139: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1140: @*/
1141: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1142: {
1146: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1147: PetscSectionGetOffset(s->field[field], point, offset);
1148: return(0);
1149: }
1153: /*@
1154: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1156: Not collective
1158: Input Parameters:
1159: + s - the PetscSection
1160: . point - the point
1161: . field - the field
1162: - offset - the offset
1164: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1166: Level: intermediate
1168: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1169: @*/
1170: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1171: {
1175: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1176: PetscSectionSetOffset(s->field[field], point, offset);
1177: return(0);
1178: }
1182: /* This gives the offset on a point of the field, ignoring constraints */
1183: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1184: {
1185: PetscInt off, foff;
1189: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1190: PetscSectionGetOffset(s, point, &off);
1191: PetscSectionGetOffset(s->field[field], point, &foff);
1192: *offset = foff - off;
1193: return(0);
1194: }
1198: /*@
1199: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1201: Not collective
1203: Input Parameter:
1204: . s - the PetscSection
1206: Output Parameters:
1207: + start - the minimum offset
1208: - end - one more than the maximum offset
1210: Level: intermediate
1212: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1213: @*/
1214: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1215: {
1216: PetscInt os = s->atlasOff[0], oe = s->atlasOff[0], pStart, pEnd, p;
1220: PetscSectionGetChart(s, &pStart, &pEnd);
1221: for (p = 0; p < pEnd-pStart; ++p) {
1222: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1224: if (off >= 0) {
1225: os = PetscMin(os, off);
1226: oe = PetscMax(oe, off+dof);
1227: }
1228: }
1229: if (start) *start = os;
1230: if (end) *end = oe;
1231: return(0);
1232: }
1236: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, PetscInt fields[], PetscSection *subs)
1237: {
1238: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1242: if (!numFields) return(0);
1243: PetscSectionGetNumFields(s, &nF);
1244: if (numFields > nF) SETERRQ2(s->atlasLayout.comm, PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1245: PetscSectionCreate(s->atlasLayout.comm, subs);
1246: PetscSectionSetNumFields(*subs, numFields);
1247: for (f = 0; f < numFields; ++f) {
1248: const char *name = NULL;
1249: PetscInt numComp = 0;
1251: PetscSectionGetFieldName(s, fields[f], &name);
1252: PetscSectionSetFieldName(*subs, f, name);
1253: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1254: PetscSectionSetFieldComponents(*subs, f, numComp);
1255: }
1256: PetscSectionGetChart(s, &pStart, &pEnd);
1257: PetscSectionSetChart(*subs, pStart, pEnd);
1258: for (p = pStart; p < pEnd; ++p) {
1259: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1261: for (f = 0; f < numFields; ++f) {
1262: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1263: PetscSectionSetFieldDof(*subs, p, f, fdof);
1264: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1265: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1266: dof += fdof;
1267: cdof += cfdof;
1268: }
1269: PetscSectionSetDof(*subs, p, dof);
1270: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1271: maxCdof = PetscMax(cdof, maxCdof);
1272: }
1273: PetscSectionSetUp(*subs);
1274: if (maxCdof) {
1275: PetscInt *indices;
1277: PetscMalloc(maxCdof * sizeof(PetscInt), &indices);
1278: for (p = pStart; p < pEnd; ++p) {
1279: PetscInt cdof;
1281: PetscSectionGetConstraintDof(*subs, p, &cdof);
1282: if (cdof) {
1283: const PetscInt *oldIndices = NULL;
1284: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1286: for (f = 0; f < numFields; ++f) {
1287: PetscInt oldFoff = 0, oldf;
1289: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1290: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1291: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1292: /* This can be sped up if we assume sorted fields */
1293: for (oldf = 0; oldf < fields[f]; ++oldf) {
1294: PetscInt oldfdof = 0;
1295: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1296: oldFoff += oldfdof;
1297: }
1298: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1299: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1300: numConst += cfdof;
1301: fOff += fdof;
1302: }
1303: PetscSectionSetConstraintIndices(*subs, p, indices);
1304: }
1305: }
1306: PetscFree(indices);
1307: }
1308: return(0);
1309: }
1313: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1314: {
1315: const PetscInt *points, *indices = NULL;
1316: PetscInt numFields, f, numSubpoints, pStart, pEnd, p, subp;
1320: PetscSectionGetNumFields(s, &numFields);
1321: PetscSectionCreate(s->atlasLayout.comm, subs);
1322: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1323: for (f = 0; f < numFields; ++f) {
1324: const char *name = NULL;
1325: PetscInt numComp = 0;
1327: PetscSectionGetFieldName(s, f, &name);
1328: PetscSectionSetFieldName(*subs, f, name);
1329: PetscSectionGetFieldComponents(s, f, &numComp);
1330: PetscSectionSetFieldComponents(*subs, f, numComp);
1331: }
1332: /* For right now, we do not try to squeeze the subchart */
1333: ISGetSize(subpointMap, &numSubpoints);
1334: ISGetIndices(subpointMap, &points);
1335: PetscSectionGetChart(s, &pStart, &pEnd);
1336: PetscSectionSetChart(*subs, 0, numSubpoints);
1337: for (p = pStart; p < pEnd; ++p) {
1338: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1340: PetscFindInt(p, numSubpoints, points, &subp);
1341: if (subp < 0) continue;
1342: for (f = 0; f < numFields; ++f) {
1343: PetscSectionGetFieldDof(s, p, f, &fdof);
1344: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1345: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1346: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1347: }
1348: PetscSectionGetDof(s, p, &dof);
1349: PetscSectionSetDof(*subs, subp, dof);
1350: PetscSectionGetConstraintDof(s, p, &cdof);
1351: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1352: }
1353: PetscSectionSetUp(*subs);
1354: /* Change offsets to original offsets */
1355: for (p = pStart; p < pEnd; ++p) {
1356: PetscInt off, foff = 0;
1358: PetscFindInt(p, numSubpoints, points, &subp);
1359: if (subp < 0) continue;
1360: for (f = 0; f < numFields; ++f) {
1361: PetscSectionGetFieldOffset(s, p, f, &foff);
1362: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1363: }
1364: PetscSectionGetOffset(s, p, &off);
1365: PetscSectionSetOffset(*subs, subp, off);
1366: }
1367: /* Copy constraint indices */
1368: for (subp = 0; subp < numSubpoints; ++subp) {
1369: PetscInt cdof;
1371: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1372: if (cdof) {
1373: for (f = 0; f < numFields; ++f) {
1374: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1375: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1376: }
1377: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1378: PetscSectionSetConstraintIndices(*subs, subp, indices);
1379: }
1380: }
1381: ISRestoreIndices(subpointMap, &points);
1382: return(0);
1383: }
1387: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1388: {
1389: PetscInt p;
1390: PetscMPIInt rank;
1394: if (s->atlasLayout.numDof != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle %d dof in a uniform section", s->atlasLayout.numDof);
1395: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1396: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
1397: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1398: for (p = 0; p < s->atlasLayout.pEnd - s->atlasLayout.pStart; ++p) {
1399: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1400: PetscInt b;
1402: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d constrained", p+s->atlasLayout.pStart, s->atlasDof[p], s->atlasOff[p]);
1403: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1404: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1405: }
1406: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1407: } else {
1408: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d\n", p+s->atlasLayout.pStart, s->atlasDof[p], s->atlasOff[p]);
1409: }
1410: }
1411: PetscViewerFlush(viewer);
1412: return(0);
1413: }
1417: /*@C
1418: PetscSectionView - Views a PetscSection
1420: Collective on PetscSection
1422: Input Parameters:
1423: + s - the PetscSection object to view
1424: - v - the viewer
1426: Level: developer
1428: .seealso PetscSectionCreate(), PetscSectionDestroy()
1429: @*/
1430: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1431: {
1432: PetscBool isascii;
1433: PetscInt f;
1437: if (!viewer) {PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer);}
1439: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1440: if (isascii) {
1441: if (s->numFields) {
1442: PetscViewerASCIIPrintf(viewer, "PetscSection with %d fields\n", s->numFields);
1443: for (f = 0; f < s->numFields; ++f) {
1444: PetscViewerASCIIPrintf(viewer, " field %d with %d components\n", f, s->numFieldComponents[f]);
1445: PetscSectionView_ASCII(s->field[f], viewer);
1446: }
1447: } else {
1448: PetscViewerASCIIPrintf(viewer, "PetscSection\n");
1449: PetscSectionView_ASCII(s, viewer);
1450: }
1451: }
1452: return(0);
1453: }
1457: /*@
1458: PetscSectionReset - Frees all section data.
1460: Not collective
1462: Input Parameters:
1463: . s - the PetscSection
1465: Level: developer
1467: .seealso: PetscSection, PetscSectionCreate()
1468: @*/
1469: PetscErrorCode PetscSectionReset(PetscSection s)
1470: {
1471: PetscInt f;
1475: PetscFree(s->numFieldComponents);
1476: for (f = 0; f < s->numFields; ++f) {
1477: PetscSectionDestroy(&s->field[f]);
1478: PetscFree(s->fieldNames[f]);
1479: }
1480: PetscFree(s->fieldNames);
1481: PetscFree(s->field);
1482: PetscSectionDestroy(&s->bc);
1483: PetscFree(s->bcIndices);
1484: PetscFree2(s->atlasDof, s->atlasOff);
1486: s->atlasLayout.pStart = -1;
1487: s->atlasLayout.pEnd = -1;
1488: s->atlasLayout.numDof = 1;
1489: s->maxDof = 0;
1490: s->setup = PETSC_FALSE;
1491: s->numFields = 0;
1492: return(0);
1493: }
1497: /*@
1498: PetscSectionDestroy - Frees a section object and frees its range if that exists.
1500: Not collective
1502: Input Parameters:
1503: . s - the PetscSection
1505: Level: developer
1507: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1508: recommended they not be used in user codes unless you really gain something in their use.
1510: .seealso: PetscSection, PetscSectionCreate()
1511: @*/
1512: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1513: {
1517: if (!*s) return(0);
1519: if (--((PetscObject)(*s))->refct > 0) {
1520: *s = NULL;
1521: return(0);
1522: }
1523: PetscSectionReset(*s);
1524: PetscHeaderDestroy(s);
1525: return(0);
1526: }
1530: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1531: {
1532: const PetscInt p = point - s->atlasLayout.pStart;
1535: *values = &baseArray[s->atlasOff[p]];
1536: return(0);
1537: }
1541: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1542: {
1543: PetscInt *array;
1544: const PetscInt p = point - s->atlasLayout.pStart;
1545: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1546: PetscInt cDim = 0;
1550: PetscSectionGetConstraintDof(s, p, &cDim);
1551: array = &baseArray[s->atlasOff[p]];
1552: if (!cDim) {
1553: if (orientation >= 0) {
1554: const PetscInt dim = s->atlasDof[p];
1555: PetscInt i;
1557: if (mode == INSERT_VALUES) {
1558: for (i = 0; i < dim; ++i) array[i] = values[i];
1559: } else {
1560: for (i = 0; i < dim; ++i) array[i] += values[i];
1561: }
1562: } else {
1563: PetscInt offset = 0;
1564: PetscInt j = -1, field, i;
1566: for (field = 0; field < s->numFields; ++field) {
1567: const PetscInt dim = s->field[field]->atlasDof[p];
1569: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1570: offset += dim;
1571: }
1572: }
1573: } else {
1574: if (orientation >= 0) {
1575: const PetscInt dim = s->atlasDof[p];
1576: PetscInt cInd = 0, i;
1577: const PetscInt *cDof;
1579: PetscSectionGetConstraintIndices(s, point, &cDof);
1580: if (mode == INSERT_VALUES) {
1581: for (i = 0; i < dim; ++i) {
1582: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1583: array[i] = values[i];
1584: }
1585: } else {
1586: for (i = 0; i < dim; ++i) {
1587: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1588: array[i] += values[i];
1589: }
1590: }
1591: } else {
1592: const PetscInt *cDof;
1593: PetscInt offset = 0;
1594: PetscInt cOffset = 0;
1595: PetscInt j = 0, field;
1597: PetscSectionGetConstraintIndices(s, point, &cDof);
1598: for (field = 0; field < s->numFields; ++field) {
1599: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
1600: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1601: const PetscInt sDim = dim - tDim;
1602: PetscInt cInd = 0, i,k;
1604: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1605: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1606: array[j] = values[k];
1607: }
1608: offset += dim;
1609: cOffset += dim - tDim;
1610: }
1611: }
1612: }
1613: return(0);
1614: }
1618: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1619: {
1623: if (s->bc) {
1624: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1625: } else *indices = NULL;
1626: return(0);
1627: }
1631: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1632: {
1636: if (s->bc) {
1637: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1638: }
1639: return(0);
1640: }
1644: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1645: {
1649: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1650: PetscSectionGetConstraintIndices(s->field[field], point, indices);
1651: return(0);
1652: }
1656: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
1657: {
1661: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1662: PetscSectionSetConstraintIndices(s->field[field], point, indices);
1663: return(0);
1664: }
1666: /*
1667: I need extract and merge routines for section based on fields. This should be trivial except for updating the
1668: constraint indices.
1670: Then I need a new interface for DMCreateDecomposition that takes groups of fields and returns a real DMPlex
1671: that shares the mesh parts, and has the extracted section
1672: */
1676: PetscErrorCode PetscSFConvertPartition(PetscSF sfPart, PetscSection partSection, IS partition, ISLocalToGlobalMapping *renumbering, PetscSF *sf)
1677: {
1678: MPI_Comm comm;
1679: PetscSF sfPoints;
1680: PetscInt *partSizes,*partOffsets,p,i,numParts,numMyPoints,numPoints,count;
1681: const PetscInt *partArray;
1682: PetscSFNode *sendPoints;
1683: PetscMPIInt rank;
1687: PetscObjectGetComm((PetscObject)sfPart,&comm);
1688: MPI_Comm_rank(comm, &rank);
1690: /* Get the number of parts and sizes that I have to distribute */
1691: PetscSFGetGraph(sfPart,NULL,&numParts,NULL,NULL);
1692: PetscMalloc2(numParts,PetscInt,&partSizes,numParts,PetscInt,&partOffsets);
1693: for (p=0,numPoints=0; p<numParts; p++) {
1694: PetscSectionGetDof(partSection, p, &partSizes[p]);
1695: numPoints += partSizes[p];
1696: }
1697: numMyPoints = 0;
1698: PetscSFFetchAndOpBegin(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1699: PetscSFFetchAndOpEnd(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1700: /* I will receive numMyPoints. I will send a total of numPoints, to be placed on remote procs at partOffsets */
1702: /* Create SF mapping locations (addressed through partition, as indexed leaves) to new owners (roots) */
1703: PetscMalloc(numPoints*sizeof(PetscSFNode),&sendPoints);
1704: for (p=0,count=0; p<numParts; p++) {
1705: for (i=0; i<partSizes[p]; i++) {
1706: sendPoints[count].rank = p;
1707: sendPoints[count].index = partOffsets[p]+i;
1708: count++;
1709: }
1710: }
1711: if (count != numPoints) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count %D should equal numPoints=%D",count,numPoints);
1712: PetscFree2(partSizes,partOffsets);
1713: ISGetIndices(partition,&partArray);
1714: PetscSFCreate(comm,&sfPoints);
1715: PetscSFSetFromOptions(sfPoints);
1716: PetscSFSetGraph(sfPoints,numMyPoints,numPoints,partArray,PETSC_USE_POINTER,sendPoints,PETSC_OWN_POINTER);
1718: /* Invert SF so that the new owners are leaves and the locations indexed through partition are the roots */
1719: PetscSFCreateInverseSF(sfPoints,sf);
1720: PetscSFDestroy(&sfPoints);
1721: ISRestoreIndices(partition,&partArray);
1723: /* Create the new local-to-global mapping */
1724: ISLocalToGlobalMappingCreateSF(*sf,0,renumbering);
1725: return(0);
1726: }
1730: /*@C
1731: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
1733: Collective
1735: Input Parameters:
1736: + sf - The SF
1737: - rootSection - Section defined on root space
1739: Output Parameters:
1740: + remoteOffsets - root offsets in leaf storage, or NULL
1741: - leafSection - Section defined on the leaf space
1743: Level: intermediate
1745: .seealso: PetscSFCreate()
1746: @*/
1747: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1748: {
1749: PetscSF embedSF;
1750: const PetscInt *ilocal, *indices;
1751: IS selected;
1752: PetscInt numFields, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;
1756: PetscSectionGetNumFields(rootSection, &numFields);
1757: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
1758: for (f = 0; f < numFields; ++f) {
1759: PetscInt numComp = 0;
1760: PetscSectionGetFieldComponents(rootSection, f, &numComp);
1761: PetscSectionSetFieldComponents(leafSection, f, numComp);
1762: }
1763: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1764: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1765: ISGetIndices(selected, &indices);
1766: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1767: ISRestoreIndices(selected, &indices);
1768: ISDestroy(&selected);
1769: PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
1770: if (ilocal) {
1771: for (i = 0; i < nleaves; ++i) {
1772: lpStart = PetscMin(lpStart, ilocal[i]);
1773: lpEnd = PetscMax(lpEnd, ilocal[i]);
1774: }
1775: ++lpEnd;
1776: } else {
1777: lpStart = 0;
1778: lpEnd = nleaves;
1779: }
1780: PetscSectionSetChart(leafSection, lpStart, lpEnd);
1781: /* Could fuse these at the cost of a copy and extra allocation */
1782: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1783: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1784: for (f = 0; f < numFields; ++f) {
1785: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1786: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1787: }
1788: if (remoteOffsets) {
1789: PetscMalloc((lpEnd - lpStart) * sizeof(PetscInt), remoteOffsets);
1790: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1791: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1792: }
1793: PetscSFDestroy(&embedSF);
1794: PetscSectionSetUp(leafSection);
1795: return(0);
1796: }
1800: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
1801: {
1802: PetscSF embedSF;
1803: const PetscInt *indices;
1804: IS selected;
1805: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0, isSize;
1809: *remoteOffsets = NULL;
1810: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
1811: if (numRoots < 0) return(0);
1812: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1813: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1814: PetscMalloc((lpEnd - lpStart) * sizeof(PetscInt), remoteOffsets);
1815: isSize = PetscMin(numRoots, rpEnd - rpStart);
1816: ISCreateStride(PETSC_COMM_SELF, isSize, rpStart, 1, &selected);
1817: ISGetIndices(selected, &indices);
1818: #if 0
1819: PetscSFCreateEmbeddedSF(sf, isSize, indices, &embedSF);
1820: #else
1821: embedSF = sf;
1822: PetscObjectReference((PetscObject) embedSF);
1823: #endif
1824: ISRestoreIndices(selected, &indices);
1825: ISDestroy(&selected);
1826: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1827: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1828: PetscSFDestroy(&embedSF);
1829: return(0);
1830: }
1834: /*@C
1835: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
1837: Input Parameters:
1838: + sf - The SF
1839: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section), or NULL
1840: - remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
1842: Output Parameters:
1843: + leafSection - Data layout of local points for incoming data (this is the distributed section)
1844: - sectionSF - The new SF
1846: Note: Either rootSection or remoteOffsets can be specified
1848: Level: intermediate
1850: .seealso: PetscSFCreate()
1851: @*/
1852: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1853: {
1854: MPI_Comm comm;
1855: const PetscInt *localPoints;
1856: const PetscSFNode *remotePoints;
1857: PetscInt lpStart, lpEnd;
1858: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
1859: PetscInt *localIndices;
1860: PetscSFNode *remoteIndices;
1861: PetscInt i, ind;
1862: PetscErrorCode ierr;
1870: PetscObjectGetComm((PetscObject)sf,&comm);
1871: PetscSFCreate(comm, sectionSF);
1872: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1873: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
1874: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
1875: if (numRoots < 0) return(0);
1876: for (i = 0; i < numPoints; ++i) {
1877: PetscInt localPoint = localPoints ? localPoints[i] : i;
1878: PetscInt dof;
1880: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
1881: PetscSectionGetDof(leafSection, localPoint, &dof);
1882: numIndices += dof;
1883: }
1884: }
1885: PetscMalloc(numIndices * sizeof(PetscInt), &localIndices);
1886: PetscMalloc(numIndices * sizeof(PetscSFNode), &remoteIndices);
1887: /* Create new index graph */
1888: for (i = 0, ind = 0; i < numPoints; ++i) {
1889: PetscInt localPoint = localPoints ? localPoints[i] : i;
1890: PetscInt rank = remotePoints[i].rank;
1892: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
1893: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
1894: PetscInt loff, dof, d;
1896: PetscSectionGetOffset(leafSection, localPoint, &loff);
1897: PetscSectionGetDof(leafSection, localPoint, &dof);
1898: for (d = 0; d < dof; ++d, ++ind) {
1899: localIndices[ind] = loff+d;
1900: remoteIndices[ind].rank = rank;
1901: remoteIndices[ind].index = remoteOffset+d;
1902: }
1903: }
1904: }
1905: PetscFree(remoteOffsets);
1906: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
1907: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
1908: return(0);
1909: }