Actual source code: vsectionis.c
petsc-3.7.7 2017-09-25
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: $ PetscSectionSetNumFields(PetscSection, numFields);
27: $ PetscSectionSetChart(PetscSection,low,high);
28: $ PetscSectionSetDof(PetscSection,point,numdof);
29: $ PetscSectionSetUp(PetscSection);
30: $ PetscSectionGetOffset(PetscSection,point,PetscInt *);
31: $ PetscSectionDestroy(PetscSection);
33: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
34: recommended they not be used in user codes unless you really gain something in their use.
36: .seealso: PetscSection, PetscSectionDestroy()
37: @*/
38: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
39: {
44: ISInitializePackage();
46: PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
48: (*s)->pStart = -1;
49: (*s)->pEnd = -1;
50: (*s)->perm = NULL;
51: (*s)->maxDof = 0;
52: (*s)->atlasDof = NULL;
53: (*s)->atlasOff = NULL;
54: (*s)->bc = NULL;
55: (*s)->bcIndices = NULL;
56: (*s)->setup = PETSC_FALSE;
57: (*s)->numFields = 0;
58: (*s)->fieldNames = NULL;
59: (*s)->field = NULL;
60: (*s)->clObj = NULL;
61: (*s)->clSection = NULL;
62: (*s)->clPoints = NULL;
63: return(0);
64: }
68: /*@
69: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
71: Collective on MPI_Comm
73: Input Parameter:
74: . section - the PetscSection
76: Output Parameter:
77: . newSection - the copy
79: Level: developer
81: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
82: @*/
83: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
84: {
85: IS perm;
86: PetscInt numFields, f, pStart, pEnd, p;
90: PetscSectionGetNumFields(section, &numFields);
91: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
92: for (f = 0; f < numFields; ++f) {
93: const char *name = NULL;
94: PetscInt numComp = 0;
96: PetscSectionGetFieldName(section, f, &name);
97: PetscSectionSetFieldName(newSection, f, name);
98: PetscSectionGetFieldComponents(section, f, &numComp);
99: PetscSectionSetFieldComponents(newSection, f, numComp);
100: }
101: PetscSectionGetChart(section, &pStart, &pEnd);
102: PetscSectionSetChart(newSection, pStart, pEnd);
103: PetscSectionGetPermutation(section, &perm);
104: PetscSectionSetPermutation(newSection, perm);
105: for (p = pStart; p < pEnd; ++p) {
106: PetscInt dof, cdof, fcdof = 0;
108: PetscSectionGetDof(section, p, &dof);
109: PetscSectionSetDof(newSection, p, dof);
110: PetscSectionGetConstraintDof(section, p, &cdof);
111: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
112: for (f = 0; f < numFields; ++f) {
113: PetscSectionGetFieldDof(section, p, f, &dof);
114: PetscSectionSetFieldDof(newSection, p, f, dof);
115: if (cdof) {
116: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
117: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
118: }
119: }
120: }
121: PetscSectionSetUp(newSection);
122: for (p = pStart; p < pEnd; ++p) {
123: PetscInt off, cdof, fcdof = 0;
124: const PetscInt *cInd;
126: /* Must set offsets in case they do not agree with the prefix sums */
127: PetscSectionGetOffset(section, p, &off);
128: PetscSectionSetOffset(newSection, p, off);
129: PetscSectionGetConstraintDof(section, p, &cdof);
130: if (cdof) {
131: PetscSectionGetConstraintIndices(section, p, &cInd);
132: PetscSectionSetConstraintIndices(newSection, p, cInd);
133: for (f = 0; f < numFields; ++f) {
134: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
135: if (fcdof) {
136: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
137: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
138: }
139: }
140: }
141: }
142: return(0);
143: }
147: /*@
148: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
150: Collective on MPI_Comm
152: Input Parameter:
153: . section - the PetscSection
155: Output Parameter:
156: . newSection - the copy
158: Level: developer
160: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
161: @*/
162: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
163: {
167: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
168: PetscSectionCopy(section, *newSection);
169: return(0);
170: }
174: /*@
175: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
177: Not collective
179: Input Parameter:
180: . s - the PetscSection
182: Output Parameter:
183: . numFields - the number of fields defined, or 0 if none were defined
185: Level: intermediate
187: .seealso: PetscSectionSetNumFields()
188: @*/
189: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
190: {
193: *numFields = s->numFields;
194: return(0);
195: }
199: /*@
200: PetscSectionSetNumFields - Sets the number of fields.
202: Not collective
204: Input Parameters:
205: + s - the PetscSection
206: - numFields - the number of fields
208: Level: intermediate
210: .seealso: PetscSectionGetNumFields()
211: @*/
212: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
213: {
214: PetscInt f;
218: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %d must be positive", numFields);
219: PetscSectionReset(s);
221: s->numFields = numFields;
222: PetscMalloc1(s->numFields, &s->numFieldComponents);
223: PetscMalloc1(s->numFields, &s->fieldNames);
224: PetscMalloc1(s->numFields, &s->field);
225: for (f = 0; f < s->numFields; ++f) {
226: char name[64];
228: s->numFieldComponents[f] = 1;
230: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
231: PetscSNPrintf(name, 64, "Field_%D", f);
232: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
233: }
234: return(0);
235: }
239: /*@C
240: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
242: Not Collective
244: Input Parameters:
245: + s - the PetscSection
246: - field - the field number
248: Output Parameter:
249: . fieldName - the field name
251: Level: developer
253: .seealso: PetscSectionSetFieldName()
254: @*/
255: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
256: {
259: 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);
260: *fieldName = s->fieldNames[field];
261: return(0);
262: }
266: /*@C
267: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
269: Not Collective
271: Input Parameters:
272: + s - the PetscSection
273: . field - the field number
274: - fieldName - the field name
276: Level: developer
278: .seealso: PetscSectionGetFieldName()
279: @*/
280: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
281: {
286: 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);
287: PetscFree(s->fieldNames[field]);
288: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
289: return(0);
290: }
294: /*@
295: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
297: Not collective
299: Input Parameters:
300: + s - the PetscSection
301: - field - the field number
303: Output Parameter:
304: . numComp - the number of field components
306: Level: intermediate
308: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
309: @*/
310: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
311: {
314: 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);
315: *numComp = s->numFieldComponents[field];
316: return(0);
317: }
321: /*@
322: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
324: Not collective
326: Input Parameters:
327: + s - the PetscSection
328: . field - the field number
329: - numComp - the number of field components
331: Level: intermediate
333: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
334: @*/
335: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
336: {
338: 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);
339: s->numFieldComponents[field] = numComp;
340: return(0);
341: }
345: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
346: {
350: if (!s->bc) {
351: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
352: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
353: }
354: return(0);
355: }
359: /*@
360: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points in the lie.
362: Not collective
364: Input Parameter:
365: . s - the PetscSection
367: Output Parameters:
368: + pStart - the first point
369: - pEnd - one past the last point
371: Level: intermediate
373: .seealso: PetscSectionSetChart(), PetscSectionCreate()
374: @*/
375: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
376: {
379: if (pStart) *pStart = s->pStart;
380: if (pEnd) *pEnd = s->pEnd;
381: return(0);
382: }
386: /*@
387: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points in the lie.
389: Not collective
391: Input Parameters:
392: + s - the PetscSection
393: . pStart - the first point
394: - pEnd - one past the last point
396: Level: intermediate
398: .seealso: PetscSectionGetChart(), PetscSectionCreate()
399: @*/
400: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
401: {
402: PetscInt f;
407: /* Cannot Reset() because it destroys field information */
408: s->setup = PETSC_FALSE;
409: PetscSectionDestroy(&s->bc);
410: PetscFree(s->bcIndices);
411: PetscFree2(s->atlasDof, s->atlasOff);
413: s->pStart = pStart;
414: s->pEnd = pEnd;
415: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
416: PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
417: for (f = 0; f < s->numFields; ++f) {
418: PetscSectionSetChart(s->field[f], pStart, pEnd);
419: }
420: return(0);
421: }
425: /*@
426: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
428: Not collective
430: Input Parameter:
431: . s - the PetscSection
433: Output Parameters:
434: . perm - The permutation as an IS
436: Level: intermediate
438: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
439: @*/
440: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
441: {
445: return(0);
446: }
450: /*@
451: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
453: Not collective
455: Input Parameters:
456: + s - the PetscSection
457: - perm - the permutation of points
459: Level: intermediate
461: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
462: @*/
463: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
464: {
468: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
469: if (s->perm != perm) {
470: ISDestroy(&s->perm);
471: s->perm = perm;
472: PetscObjectReference((PetscObject) s->perm);
473: }
474: return(0);
475: }
479: /*@
480: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
482: Not collective
484: Input Parameters:
485: + s - the PetscSection
486: - point - the point
488: Output Parameter:
489: . numDof - the number of dof
491: Level: intermediate
493: .seealso: PetscSectionSetDof(), PetscSectionCreate()
494: @*/
495: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
496: {
498: #if defined(PETSC_USE_DEBUG)
499: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
500: #endif
501: *numDof = s->atlasDof[point - s->pStart];
502: return(0);
503: }
507: /*@
508: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
510: Not collective
512: Input Parameters:
513: + s - the PetscSection
514: . point - the point
515: - numDof - the number of dof
517: Level: intermediate
519: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
520: @*/
521: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
522: {
524: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
525: s->atlasDof[point - s->pStart] = numDof;
526: return(0);
527: }
531: /*@
532: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
534: Not collective
536: Input Parameters:
537: + s - the PetscSection
538: . point - the point
539: - numDof - the number of additional dof
541: Level: intermediate
543: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
544: @*/
545: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
546: {
548: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
549: s->atlasDof[point - s->pStart] += numDof;
550: return(0);
551: }
555: /*@
556: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
558: Not collective
560: Input Parameters:
561: + s - the PetscSection
562: . point - the point
563: - field - the field
565: Output Parameter:
566: . numDof - the number of dof
568: Level: intermediate
570: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
571: @*/
572: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
573: {
577: 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);
578: PetscSectionGetDof(s->field[field], point, numDof);
579: return(0);
580: }
584: /*@
585: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
587: Not collective
589: Input Parameters:
590: + s - the PetscSection
591: . point - the point
592: . field - the field
593: - numDof - the number of dof
595: Level: intermediate
597: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
598: @*/
599: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
600: {
604: 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);
605: PetscSectionSetDof(s->field[field], point, numDof);
606: return(0);
607: }
611: /*@
612: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
614: Not collective
616: Input Parameters:
617: + s - the PetscSection
618: . point - the point
619: . field - the field
620: - numDof - the number of dof
622: Level: intermediate
624: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
625: @*/
626: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
627: {
631: 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);
632: PetscSectionAddDof(s->field[field], point, numDof);
633: return(0);
634: }
638: /*@
639: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
641: Not collective
643: Input Parameters:
644: + s - the PetscSection
645: - point - the point
647: Output Parameter:
648: . numDof - the number of dof which are fixed by constraints
650: Level: intermediate
652: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
653: @*/
654: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
655: {
659: if (s->bc) {
660: PetscSectionGetDof(s->bc, point, numDof);
661: } else *numDof = 0;
662: return(0);
663: }
667: /*@
668: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
670: Not collective
672: Input Parameters:
673: + s - the PetscSection
674: . point - the point
675: - numDof - the number of dof which are fixed by constraints
677: Level: intermediate
679: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
680: @*/
681: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
682: {
686: if (numDof) {
687: PetscSectionCheckConstraints_Static(s);
688: PetscSectionSetDof(s->bc, point, numDof);
689: }
690: return(0);
691: }
695: /*@
696: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
698: Not collective
700: Input Parameters:
701: + s - the PetscSection
702: . point - the point
703: - numDof - the number of additional dof which are fixed by constraints
705: Level: intermediate
707: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
708: @*/
709: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
710: {
714: if (numDof) {
715: PetscSectionCheckConstraints_Static(s);
716: PetscSectionAddDof(s->bc, point, numDof);
717: }
718: return(0);
719: }
723: /*@
724: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
726: Not collective
728: Input Parameters:
729: + s - the PetscSection
730: . point - the point
731: - field - the field
733: Output Parameter:
734: . numDof - the number of dof which are fixed by constraints
736: Level: intermediate
738: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
739: @*/
740: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
741: {
745: 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);
746: PetscSectionGetConstraintDof(s->field[field], point, numDof);
747: return(0);
748: }
752: /*@
753: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
755: Not collective
757: Input Parameters:
758: + s - the PetscSection
759: . point - the point
760: . field - the field
761: - numDof - the number of dof which are fixed by constraints
763: Level: intermediate
765: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
766: @*/
767: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
768: {
772: 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);
773: PetscSectionSetConstraintDof(s->field[field], point, numDof);
774: return(0);
775: }
779: /*@
780: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
782: Not collective
784: Input Parameters:
785: + s - the PetscSection
786: . point - the point
787: . field - the field
788: - numDof - the number of additional dof which are fixed by constraints
790: Level: intermediate
792: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
793: @*/
794: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
795: {
799: 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);
800: PetscSectionAddConstraintDof(s->field[field], point, numDof);
801: return(0);
802: }
806: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
807: {
811: if (s->bc) {
812: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
814: PetscSectionSetUp(s->bc);
815: PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
816: }
817: return(0);
818: }
822: /*@
823: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
825: Not collective
827: Input Parameter:
828: . s - the PetscSection
830: Level: intermediate
832: .seealso: PetscSectionCreate()
833: @*/
834: PetscErrorCode PetscSectionSetUp(PetscSection s)
835: {
836: const PetscInt *pind = NULL;
837: PetscInt offset = 0, p, f;
838: PetscErrorCode ierr;
841: if (s->setup) return(0);
842: s->setup = PETSC_TRUE;
843: if (s->perm) {ISGetIndices(s->perm, &pind);}
844: for (p = 0; p < s->pEnd - s->pStart; ++p) {
845: const PetscInt q = pind ? pind[p] : p;
847: s->atlasOff[q] = offset;
848: offset += s->atlasDof[q];
849: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
850: }
851: PetscSectionSetUpBC(s);
852: /* Assume that all fields have the same chart */
853: for (p = 0; p < s->pEnd - s->pStart; ++p) {
854: const PetscInt q = pind ? pind[p] : p;
855: PetscInt off = s->atlasOff[q];
857: for (f = 0; f < s->numFields; ++f) {
858: PetscSection sf = s->field[f];
860: sf->atlasOff[q] = off;
861: off += sf->atlasDof[q];
862: }
863: }
864: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
865: for (f = 0; f < s->numFields; ++f) {
866: PetscSectionSetUpBC(s->field[f]);
867: }
868: return(0);
869: }
873: /*@
874: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
876: Not collective
878: Input Parameters:
879: . s - the PetscSection
881: Output Parameter:
882: . maxDof - the maximum dof
884: Level: intermediate
886: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
887: @*/
888: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
889: {
891: *maxDof = s->maxDof;
892: return(0);
893: }
897: /*@
898: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
900: Not collective
902: Input Parameters:
903: + s - the PetscSection
904: - size - the allocated size
906: Output Parameter:
907: . size - the size of an array which can hold all the dofs
909: Level: intermediate
911: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
912: @*/
913: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
914: {
915: PetscInt p, n = 0;
918: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
919: *size = n;
920: return(0);
921: }
925: /*@
926: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
928: Not collective
930: Input Parameters:
931: + s - the PetscSection
932: - point - the point
934: Output Parameter:
935: . size - the size of an array which can hold all unconstrained dofs
937: Level: intermediate
939: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
940: @*/
941: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
942: {
943: PetscInt p, n = 0;
946: for (p = 0; p < s->pEnd - s->pStart; ++p) {
947: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
948: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
949: }
950: *size = n;
951: return(0);
952: }
956: /*@
957: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
958: the local section and an SF describing the section point overlap.
960: Input Parameters:
961: + s - The PetscSection for the local field layout
962: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
963: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
964: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
966: Output Parameter:
967: . gsection - The PetscSection for the global field layout
969: Note: This gives negative sizes and offsets to points not owned by this process
971: Level: developer
973: .seealso: PetscSectionCreate()
974: @*/
975: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
976: {
977: const PetscInt *pind = NULL;
978: PetscInt *recv = NULL, *neg = NULL;
979: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
980: PetscErrorCode ierr;
983: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
984: PetscSectionGetChart(s, &pStart, &pEnd);
985: PetscSectionSetChart(*gsection, pStart, pEnd);
986: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
987: nlocal = nroots; /* The local/leaf space matches global/root space */
988: /* Must allocate for all points visible to SF, which may be more than this section */
989: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
990: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
991: PetscMalloc2(nroots,&neg,nlocal,&recv);
992: PetscMemzero(neg,nroots*sizeof(PetscInt));
993: }
994: /* Mark all local points with negative dof */
995: for (p = pStart; p < pEnd; ++p) {
996: PetscSectionGetDof(s, p, &dof);
997: PetscSectionSetDof(*gsection, p, dof);
998: PetscSectionGetConstraintDof(s, p, &cdof);
999: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1000: if (neg) neg[p] = -(dof+1);
1001: }
1002: PetscSectionSetUpBC(*gsection);
1003: if (nroots >= 0) {
1004: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1005: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1006: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1007: for (p = pStart; p < pEnd; ++p) {
1008: if (recv[p] < 0) {
1009: (*gsection)->atlasDof[p-pStart] = recv[p];
1010: PetscSectionGetDof(s, p, &dof);
1011: if (-(recv[p]+1) != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %d for point %d is not the unconstrained %d", -(recv[p]+1), p, dof);
1012: }
1013: }
1014: }
1015: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1016: if (s->perm) {ISGetIndices(s->perm, &pind);}
1017: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1018: const PetscInt q = pind ? pind[p] : p;
1020: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1021: (*gsection)->atlasOff[q] = off;
1022: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1023: }
1024: if (!localOffsets) {
1025: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1026: globalOff -= off;
1027: }
1028: for (p = pStart, off = 0; p < pEnd; ++p) {
1029: (*gsection)->atlasOff[p-pStart] += globalOff;
1030: if (neg) neg[p] = -((*gsection)->atlasOff[p-pStart]+1);
1031: }
1032: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1033: /* Put in negative offsets for ghost points */
1034: if (nroots >= 0) {
1035: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1036: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1037: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1038: for (p = pStart; p < pEnd; ++p) {
1039: if (recv[p] < 0) (*gsection)->atlasOff[p-pStart] = recv[p];
1040: }
1041: }
1042: PetscFree2(neg,recv);
1043: PetscSectionViewFromOptions(*gsection,NULL,"-global_section_view");
1044: return(0);
1045: }
1049: /*@
1050: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1051: the local section and an SF describing the section point overlap.
1053: Input Parameters:
1054: + s - The PetscSection for the local field layout
1055: . sf - The SF describing parallel layout of the section points
1056: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1057: . numExcludes - The number of exclusion ranges
1058: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1060: Output Parameter:
1061: . gsection - The PetscSection for the global field layout
1063: Note: This gives negative sizes and offsets to points not owned by this process
1065: Level: developer
1067: .seealso: PetscSectionCreate()
1068: @*/
1069: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1070: {
1071: const PetscInt *pind = NULL;
1072: PetscInt *neg = NULL, *tmpOff = NULL;
1073: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1074: PetscErrorCode ierr;
1077: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1078: PetscSectionGetChart(s, &pStart, &pEnd);
1079: PetscSectionSetChart(*gsection, pStart, pEnd);
1080: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1081: if (nroots >= 0) {
1082: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1083: PetscCalloc1(nroots, &neg);
1084: if (nroots > pEnd-pStart) {
1085: PetscCalloc1(nroots, &tmpOff);
1086: } else {
1087: tmpOff = &(*gsection)->atlasDof[-pStart];
1088: }
1089: }
1090: /* Mark ghost points with negative dof */
1091: for (p = pStart; p < pEnd; ++p) {
1092: for (e = 0; e < numExcludes; ++e) {
1093: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1094: PetscSectionSetDof(*gsection, p, 0);
1095: break;
1096: }
1097: }
1098: if (e < numExcludes) continue;
1099: PetscSectionGetDof(s, p, &dof);
1100: PetscSectionSetDof(*gsection, p, dof);
1101: PetscSectionGetConstraintDof(s, p, &cdof);
1102: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1103: if (neg) neg[p] = -(dof+1);
1104: }
1105: PetscSectionSetUpBC(*gsection);
1106: if (nroots >= 0) {
1107: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1108: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1109: if (nroots > pEnd - pStart) {
1110: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1111: }
1112: }
1113: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1114: if (s->perm) {ISGetIndices(s->perm, &pind);}
1115: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1116: const PetscInt q = pind ? pind[p] : p;
1118: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1119: (*gsection)->atlasOff[q] = off;
1120: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1121: }
1122: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1123: globalOff -= off;
1124: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1125: (*gsection)->atlasOff[p] += globalOff;
1126: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1127: }
1128: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1129: /* Put in negative offsets for ghost points */
1130: if (nroots >= 0) {
1131: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1132: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1133: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1134: if (nroots > pEnd - pStart) {
1135: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1136: }
1137: }
1138: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1139: PetscFree(neg);
1140: return(0);
1141: }
1145: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1146: {
1147: PetscInt pStart, pEnd, p, localSize = 0;
1151: PetscSectionGetChart(s, &pStart, &pEnd);
1152: for (p = pStart; p < pEnd; ++p) {
1153: PetscInt dof;
1155: PetscSectionGetDof(s, p, &dof);
1156: if (dof > 0) ++localSize;
1157: }
1158: PetscLayoutCreate(comm, layout);
1159: PetscLayoutSetLocalSize(*layout, localSize);
1160: PetscLayoutSetBlockSize(*layout, 1);
1161: PetscLayoutSetUp(*layout);
1162: return(0);
1163: }
1167: /*@
1168: PetscSectionGetValueLayout - Get the PetscLayout associated with a section, usually the default global section.
1170: Input Parameters:
1171: + comm - The MPI_Comm
1172: - s - The PetscSection
1174: Output Parameter:
1175: . layout - The layout for the section
1177: Level: developer
1179: .seealso: PetscSectionCreate()
1180: @*/
1181: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1182: {
1183: PetscInt pStart, pEnd, p, localSize = 0;
1187: PetscSectionGetChart(s, &pStart, &pEnd);
1188: for (p = pStart; p < pEnd; ++p) {
1189: PetscInt dof,cdof;
1191: PetscSectionGetDof(s, p, &dof);
1192: PetscSectionGetConstraintDof(s, p, &cdof);
1193: if (dof-cdof > 0) localSize += dof-cdof;
1194: }
1195: PetscLayoutCreate(comm, layout);
1196: PetscLayoutSetLocalSize(*layout, localSize);
1197: PetscLayoutSetBlockSize(*layout, 1);
1198: PetscLayoutSetUp(*layout);
1199: return(0);
1200: }
1204: /*@
1205: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1207: Not collective
1209: Input Parameters:
1210: + s - the PetscSection
1211: - point - the point
1213: Output Parameter:
1214: . offset - the offset
1216: Level: intermediate
1218: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1219: @*/
1220: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1221: {
1223: #if defined(PETSC_USE_DEBUG)
1224: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
1225: #endif
1226: *offset = s->atlasOff[point - s->pStart];
1227: return(0);
1228: }
1232: /*@
1233: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1235: Not collective
1237: Input Parameters:
1238: + s - the PetscSection
1239: . point - the point
1240: - offset - the offset
1242: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1244: Level: intermediate
1246: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1247: @*/
1248: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1249: {
1251: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
1252: s->atlasOff[point - s->pStart] = offset;
1253: return(0);
1254: }
1258: /*@
1259: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1261: Not collective
1263: Input Parameters:
1264: + s - the PetscSection
1265: . point - the point
1266: - field - the field
1268: Output Parameter:
1269: . offset - the offset
1271: Level: intermediate
1273: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1274: @*/
1275: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1276: {
1280: 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);
1281: PetscSectionGetOffset(s->field[field], point, offset);
1282: return(0);
1283: }
1287: /*@
1288: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1290: Not collective
1292: Input Parameters:
1293: + s - the PetscSection
1294: . point - the point
1295: . field - the field
1296: - offset - the offset
1298: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1300: Level: intermediate
1302: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1303: @*/
1304: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1305: {
1309: 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);
1310: PetscSectionSetOffset(s->field[field], point, offset);
1311: return(0);
1312: }
1316: /* This gives the offset on a point of the field, ignoring constraints */
1317: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1318: {
1319: PetscInt off, foff;
1323: 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);
1324: PetscSectionGetOffset(s, point, &off);
1325: PetscSectionGetOffset(s->field[field], point, &foff);
1326: *offset = foff - off;
1327: return(0);
1328: }
1332: /*@
1333: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1335: Not collective
1337: Input Parameter:
1338: . s - the PetscSection
1340: Output Parameters:
1341: + start - the minimum offset
1342: - end - one more than the maximum offset
1344: Level: intermediate
1346: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1347: @*/
1348: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1349: {
1350: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1354: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1355: PetscSectionGetChart(s, &pStart, &pEnd);
1356: for (p = 0; p < pEnd-pStart; ++p) {
1357: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1359: if (off >= 0) {
1360: os = PetscMin(os, off);
1361: oe = PetscMax(oe, off+dof);
1362: }
1363: }
1364: if (start) *start = os;
1365: if (end) *end = oe;
1366: return(0);
1367: }
1371: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, PetscInt fields[], PetscSection *subs)
1372: {
1373: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1377: if (!numFields) return(0);
1378: PetscSectionGetNumFields(s, &nF);
1379: if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1380: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1381: PetscSectionSetNumFields(*subs, numFields);
1382: for (f = 0; f < numFields; ++f) {
1383: const char *name = NULL;
1384: PetscInt numComp = 0;
1386: PetscSectionGetFieldName(s, fields[f], &name);
1387: PetscSectionSetFieldName(*subs, f, name);
1388: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1389: PetscSectionSetFieldComponents(*subs, f, numComp);
1390: }
1391: PetscSectionGetChart(s, &pStart, &pEnd);
1392: PetscSectionSetChart(*subs, pStart, pEnd);
1393: for (p = pStart; p < pEnd; ++p) {
1394: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1396: for (f = 0; f < numFields; ++f) {
1397: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1398: PetscSectionSetFieldDof(*subs, p, f, fdof);
1399: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1400: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1401: dof += fdof;
1402: cdof += cfdof;
1403: }
1404: PetscSectionSetDof(*subs, p, dof);
1405: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1406: maxCdof = PetscMax(cdof, maxCdof);
1407: }
1408: PetscSectionSetUp(*subs);
1409: if (maxCdof) {
1410: PetscInt *indices;
1412: PetscMalloc1(maxCdof, &indices);
1413: for (p = pStart; p < pEnd; ++p) {
1414: PetscInt cdof;
1416: PetscSectionGetConstraintDof(*subs, p, &cdof);
1417: if (cdof) {
1418: const PetscInt *oldIndices = NULL;
1419: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1421: for (f = 0; f < numFields; ++f) {
1422: PetscInt oldFoff = 0, oldf;
1424: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1425: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1426: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1427: /* This can be sped up if we assume sorted fields */
1428: for (oldf = 0; oldf < fields[f]; ++oldf) {
1429: PetscInt oldfdof = 0;
1430: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1431: oldFoff += oldfdof;
1432: }
1433: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1434: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1435: numConst += cfdof;
1436: fOff += fdof;
1437: }
1438: PetscSectionSetConstraintIndices(*subs, p, indices);
1439: }
1440: }
1441: PetscFree(indices);
1442: }
1443: return(0);
1444: }
1448: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1449: {
1450: const PetscInt *points = NULL, *indices = NULL;
1451: PetscInt numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;
1455: PetscSectionGetNumFields(s, &numFields);
1456: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1457: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1458: for (f = 0; f < numFields; ++f) {
1459: const char *name = NULL;
1460: PetscInt numComp = 0;
1462: PetscSectionGetFieldName(s, f, &name);
1463: PetscSectionSetFieldName(*subs, f, name);
1464: PetscSectionGetFieldComponents(s, f, &numComp);
1465: PetscSectionSetFieldComponents(*subs, f, numComp);
1466: }
1467: /* For right now, we do not try to squeeze the subchart */
1468: if (subpointMap) {
1469: ISGetSize(subpointMap, &numSubpoints);
1470: ISGetIndices(subpointMap, &points);
1471: }
1472: PetscSectionGetChart(s, &pStart, &pEnd);
1473: PetscSectionSetChart(*subs, 0, numSubpoints);
1474: for (p = pStart; p < pEnd; ++p) {
1475: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1477: PetscFindInt(p, numSubpoints, points, &subp);
1478: if (subp < 0) continue;
1479: for (f = 0; f < numFields; ++f) {
1480: PetscSectionGetFieldDof(s, p, f, &fdof);
1481: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1482: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1483: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1484: }
1485: PetscSectionGetDof(s, p, &dof);
1486: PetscSectionSetDof(*subs, subp, dof);
1487: PetscSectionGetConstraintDof(s, p, &cdof);
1488: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1489: }
1490: PetscSectionSetUp(*subs);
1491: /* Change offsets to original offsets */
1492: for (p = pStart; p < pEnd; ++p) {
1493: PetscInt off, foff = 0;
1495: PetscFindInt(p, numSubpoints, points, &subp);
1496: if (subp < 0) continue;
1497: for (f = 0; f < numFields; ++f) {
1498: PetscSectionGetFieldOffset(s, p, f, &foff);
1499: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1500: }
1501: PetscSectionGetOffset(s, p, &off);
1502: PetscSectionSetOffset(*subs, subp, off);
1503: }
1504: /* Copy constraint indices */
1505: for (subp = 0; subp < numSubpoints; ++subp) {
1506: PetscInt cdof;
1508: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1509: if (cdof) {
1510: for (f = 0; f < numFields; ++f) {
1511: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1512: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1513: }
1514: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1515: PetscSectionSetConstraintIndices(*subs, subp, indices);
1516: }
1517: }
1518: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1519: return(0);
1520: }
1524: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1525: {
1526: PetscInt p;
1527: PetscMPIInt rank;
1531: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1532: PetscViewerASCIIPushSynchronized(viewer);
1533: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1534: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1535: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1536: PetscInt b;
1538: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1539: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1540: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1541: }
1542: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1543: } else {
1544: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1545: }
1546: }
1547: PetscViewerFlush(viewer);
1548: PetscViewerASCIIPopSynchronized(viewer);
1549: return(0);
1550: }
1554: /*@C
1555: PetscSectionView - Views a PetscSection
1557: Collective on PetscSection
1559: Input Parameters:
1560: + s - the PetscSection object to view
1561: - v - the viewer
1563: Level: developer
1565: .seealso PetscSectionCreate(), PetscSectionDestroy()
1566: @*/
1567: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1568: {
1569: PetscBool isascii;
1570: PetscInt f;
1575: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1577: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1578: if (isascii) {
1579: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1580: if (s->numFields) {
1581: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1582: for (f = 0; f < s->numFields; ++f) {
1583: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
1584: PetscSectionView_ASCII(s->field[f], viewer);
1585: }
1586: } else {
1587: PetscSectionView_ASCII(s, viewer);
1588: }
1589: }
1590: return(0);
1591: }
1595: /*@
1596: PetscSectionReset - Frees all section data.
1598: Not collective
1600: Input Parameters:
1601: . s - the PetscSection
1603: Level: developer
1605: .seealso: PetscSection, PetscSectionCreate()
1606: @*/
1607: PetscErrorCode PetscSectionReset(PetscSection s)
1608: {
1609: PetscInt f;
1613: PetscFree(s->numFieldComponents);
1614: for (f = 0; f < s->numFields; ++f) {
1615: PetscSectionDestroy(&s->field[f]);
1616: PetscFree(s->fieldNames[f]);
1617: }
1618: PetscFree(s->fieldNames);
1619: PetscFree(s->field);
1620: PetscSectionDestroy(&s->bc);
1621: PetscFree(s->bcIndices);
1622: PetscFree2(s->atlasDof, s->atlasOff);
1623: PetscSectionDestroy(&s->clSection);
1624: ISDestroy(&s->clPoints);
1625: ISDestroy(&s->perm);
1627: s->pStart = -1;
1628: s->pEnd = -1;
1629: s->maxDof = 0;
1630: s->setup = PETSC_FALSE;
1631: s->numFields = 0;
1632: s->clObj = NULL;
1633: return(0);
1634: }
1638: /*@
1639: PetscSectionDestroy - Frees a section object and frees its range if that exists.
1641: Not collective
1643: Input Parameters:
1644: . s - the PetscSection
1646: Level: developer
1648: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1649: recommended they not be used in user codes unless you really gain something in their use.
1651: .seealso: PetscSection, PetscSectionCreate()
1652: @*/
1653: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1654: {
1658: if (!*s) return(0);
1660: if (--((PetscObject)(*s))->refct > 0) {
1661: *s = NULL;
1662: return(0);
1663: }
1664: PetscSectionReset(*s);
1665: PetscHeaderDestroy(s);
1666: return(0);
1667: }
1671: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1672: {
1673: const PetscInt p = point - s->pStart;
1676: *values = &baseArray[s->atlasOff[p]];
1677: return(0);
1678: }
1682: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1683: {
1684: PetscInt *array;
1685: const PetscInt p = point - s->pStart;
1686: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1687: PetscInt cDim = 0;
1691: PetscSectionGetConstraintDof(s, p, &cDim);
1692: array = &baseArray[s->atlasOff[p]];
1693: if (!cDim) {
1694: if (orientation >= 0) {
1695: const PetscInt dim = s->atlasDof[p];
1696: PetscInt i;
1698: if (mode == INSERT_VALUES) {
1699: for (i = 0; i < dim; ++i) array[i] = values[i];
1700: } else {
1701: for (i = 0; i < dim; ++i) array[i] += values[i];
1702: }
1703: } else {
1704: PetscInt offset = 0;
1705: PetscInt j = -1, field, i;
1707: for (field = 0; field < s->numFields; ++field) {
1708: const PetscInt dim = s->field[field]->atlasDof[p];
1710: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1711: offset += dim;
1712: }
1713: }
1714: } else {
1715: if (orientation >= 0) {
1716: const PetscInt dim = s->atlasDof[p];
1717: PetscInt cInd = 0, i;
1718: const PetscInt *cDof;
1720: PetscSectionGetConstraintIndices(s, point, &cDof);
1721: if (mode == INSERT_VALUES) {
1722: for (i = 0; i < dim; ++i) {
1723: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1724: array[i] = values[i];
1725: }
1726: } else {
1727: for (i = 0; i < dim; ++i) {
1728: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1729: array[i] += values[i];
1730: }
1731: }
1732: } else {
1733: const PetscInt *cDof;
1734: PetscInt offset = 0;
1735: PetscInt cOffset = 0;
1736: PetscInt j = 0, field;
1738: PetscSectionGetConstraintIndices(s, point, &cDof);
1739: for (field = 0; field < s->numFields; ++field) {
1740: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
1741: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1742: const PetscInt sDim = dim - tDim;
1743: PetscInt cInd = 0, i,k;
1745: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1746: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1747: array[j] = values[k];
1748: }
1749: offset += dim;
1750: cOffset += dim - tDim;
1751: }
1752: }
1753: }
1754: return(0);
1755: }
1759: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1760: {
1764: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1765: return(0);
1766: }
1770: /*@C
1771: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
1773: Input Parameters:
1774: + s - The PetscSection
1775: - point - The point
1777: Output Parameter:
1778: . indices - The constrained dofs
1780: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
1782: Level: advanced
1784: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1785: @*/
1786: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1787: {
1791: if (s->bc) {
1792: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1793: } else *indices = NULL;
1794: return(0);
1795: }
1799: /*@C
1800: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
1802: Input Parameters:
1803: + s - The PetscSection
1804: . point - The point
1805: - indices - The constrained dofs
1807: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
1809: Level: advanced
1811: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1812: @*/
1813: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1814: {
1818: if (s->bc) {
1819: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1820: }
1821: return(0);
1822: }
1826: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1827: {
1831: 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);
1832: PetscSectionGetConstraintIndices(s->field[field], point, indices);
1833: return(0);
1834: }
1838: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
1839: {
1843: 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);
1844: PetscSectionSetConstraintIndices(s->field[field], point, indices);
1845: return(0);
1846: }
1850: /*@
1851: PetscSectionPermute - Reorder the section according to the input point permutation
1853: Collective on PetscSection
1855: Input Parameter:
1856: + section - The PetscSection object
1857: - perm - The point permutation, old point p becomes new point perm[p]
1859: Output Parameter:
1860: . sectionNew - The permuted PetscSection
1862: Level: intermediate
1864: .keywords: mesh
1865: .seealso: MatPermute()
1866: @*/
1867: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1868: {
1869: PetscSection s = section, sNew;
1870: const PetscInt *perm;
1871: PetscInt numFields, f, numPoints, pStart, pEnd, p;
1872: PetscErrorCode ierr;
1878: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1879: PetscSectionGetNumFields(s, &numFields);
1880: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1881: for (f = 0; f < numFields; ++f) {
1882: const char *name;
1883: PetscInt numComp;
1885: PetscSectionGetFieldName(s, f, &name);
1886: PetscSectionSetFieldName(sNew, f, name);
1887: PetscSectionGetFieldComponents(s, f, &numComp);
1888: PetscSectionSetFieldComponents(sNew, f, numComp);
1889: }
1890: ISGetLocalSize(permutation, &numPoints);
1891: ISGetIndices(permutation, &perm);
1892: PetscSectionGetChart(s, &pStart, &pEnd);
1893: PetscSectionSetChart(sNew, pStart, pEnd);
1894: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1895: for (p = pStart; p < pEnd; ++p) {
1896: PetscInt dof, cdof;
1898: PetscSectionGetDof(s, p, &dof);
1899: PetscSectionSetDof(sNew, perm[p], dof);
1900: PetscSectionGetConstraintDof(s, p, &cdof);
1901: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1902: for (f = 0; f < numFields; ++f) {
1903: PetscSectionGetFieldDof(s, p, f, &dof);
1904: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1905: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1906: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1907: }
1908: }
1909: PetscSectionSetUp(sNew);
1910: for (p = pStart; p < pEnd; ++p) {
1911: const PetscInt *cind;
1912: PetscInt cdof;
1914: PetscSectionGetConstraintDof(s, p, &cdof);
1915: if (cdof) {
1916: PetscSectionGetConstraintIndices(s, p, &cind);
1917: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1918: }
1919: for (f = 0; f < numFields; ++f) {
1920: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1921: if (cdof) {
1922: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1923: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1924: }
1925: }
1926: }
1927: ISRestoreIndices(permutation, &perm);
1928: *sectionNew = sNew;
1929: return(0);
1930: }
1934: /*@C
1935: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
1937: Collective
1939: Input Parameters:
1940: + sf - The SF
1941: - rootSection - Section defined on root space
1943: Output Parameters:
1944: + remoteOffsets - root offsets in leaf storage, or NULL
1945: - leafSection - Section defined on the leaf space
1947: Level: intermediate
1949: .seealso: PetscSFCreate()
1950: @*/
1951: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1952: {
1953: PetscSF embedSF;
1954: const PetscInt *ilocal, *indices;
1955: IS selected;
1956: PetscInt numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;
1960: PetscSectionGetNumFields(rootSection, &numFields);
1961: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
1962: for (f = 0; f < numFields; ++f) {
1963: PetscInt numComp = 0;
1964: PetscSectionGetFieldComponents(rootSection, f, &numComp);
1965: PetscSectionSetFieldComponents(leafSection, f, numComp);
1966: }
1967: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1968: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
1969: rpEnd = PetscMin(rpEnd,nroots);
1970: rpEnd = PetscMax(rpStart,rpEnd);
1971: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1972: ISGetIndices(selected, &indices);
1973: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1974: ISRestoreIndices(selected, &indices);
1975: ISDestroy(&selected);
1976: PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
1977: if (nleaves && ilocal) {
1978: for (i = 0; i < nleaves; ++i) {
1979: lpStart = PetscMin(lpStart, ilocal[i]);
1980: lpEnd = PetscMax(lpEnd, ilocal[i]);
1981: }
1982: ++lpEnd;
1983: } else {
1984: lpStart = 0;
1985: lpEnd = nleaves;
1986: }
1987: PetscSectionSetChart(leafSection, lpStart, lpEnd);
1988: /* Could fuse these at the cost of a copy and extra allocation */
1989: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1990: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1991: for (f = 0; f < numFields; ++f) {
1992: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1993: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1994: }
1995: if (remoteOffsets) {
1996: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
1997: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1998: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1999: }
2000: PetscSFDestroy(&embedSF);
2001: PetscSectionSetUp(leafSection);
2002: return(0);
2003: }
2007: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2008: {
2009: PetscSF embedSF;
2010: const PetscInt *indices;
2011: IS selected;
2012: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2013: PetscErrorCode ierr;
2016: *remoteOffsets = NULL;
2017: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2018: if (numRoots < 0) return(0);
2019: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2020: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2021: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2022: ISGetIndices(selected, &indices);
2023: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2024: ISRestoreIndices(selected, &indices);
2025: ISDestroy(&selected);
2026: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2027: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2028: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2029: PetscSFDestroy(&embedSF);
2030: return(0);
2031: }
2035: /*@C
2036: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2038: Input Parameters:
2039: + sf - The SF
2040: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2041: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2042: - leafSection - Data layout of local points for incoming data (this is the distributed section)
2044: Output Parameters:
2045: - sectionSF - The new SF
2047: Note: Either rootSection or remoteOffsets can be specified
2049: Level: intermediate
2051: .seealso: PetscSFCreate()
2052: @*/
2053: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2054: {
2055: MPI_Comm comm;
2056: const PetscInt *localPoints;
2057: const PetscSFNode *remotePoints;
2058: PetscInt lpStart, lpEnd;
2059: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2060: PetscInt *localIndices;
2061: PetscSFNode *remoteIndices;
2062: PetscInt i, ind;
2063: PetscErrorCode ierr;
2071: PetscObjectGetComm((PetscObject)sf,&comm);
2072: PetscSFCreate(comm, sectionSF);
2073: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2074: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2075: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2076: if (numRoots < 0) return(0);
2077: for (i = 0; i < numPoints; ++i) {
2078: PetscInt localPoint = localPoints ? localPoints[i] : i;
2079: PetscInt dof;
2081: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2082: PetscSectionGetDof(leafSection, localPoint, &dof);
2083: numIndices += dof;
2084: }
2085: }
2086: PetscMalloc1(numIndices, &localIndices);
2087: PetscMalloc1(numIndices, &remoteIndices);
2088: /* Create new index graph */
2089: for (i = 0, ind = 0; i < numPoints; ++i) {
2090: PetscInt localPoint = localPoints ? localPoints[i] : i;
2091: PetscInt rank = remotePoints[i].rank;
2093: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2094: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2095: PetscInt loff, dof, d;
2097: PetscSectionGetOffset(leafSection, localPoint, &loff);
2098: PetscSectionGetDof(leafSection, localPoint, &dof);
2099: for (d = 0; d < dof; ++d, ++ind) {
2100: localIndices[ind] = loff+d;
2101: remoteIndices[ind].rank = rank;
2102: remoteIndices[ind].index = remoteOffset+d;
2103: }
2104: }
2105: }
2106: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2107: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2108: return(0);
2109: }
2113: /*@
2114: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2116: Input Parameters:
2117: + section - The PetscSection
2118: . obj - A PetscObject which serves as the key for this index
2119: . clSection - Section giving the size of the closure of each point
2120: - clPoints - IS giving the points in each closure
2122: Note: We compress out closure points with no dofs in this section
2124: Level: intermediate
2126: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2127: @*/
2128: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2129: {
2133: section->clObj = obj;
2134: PetscSectionDestroy(§ion->clSection);
2135: ISDestroy(§ion->clPoints);
2136: section->clSection = clSection;
2137: section->clPoints = clPoints;
2138: return(0);
2139: }
2143: /*@
2144: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2146: Input Parameters:
2147: + section - The PetscSection
2148: - obj - A PetscObject which serves as the key for this index
2150: Output Parameters:
2151: + clSection - Section giving the size of the closure of each point
2152: - clPoints - IS giving the points in each closure
2154: Note: We compress out closure points with no dofs in this section
2156: Level: intermediate
2158: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2159: @*/
2160: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2161: {
2163: if (section->clObj == obj) {
2164: if (clSection) *clSection = section->clSection;
2165: if (clPoints) *clPoints = section->clPoints;
2166: } else {
2167: if (clSection) *clSection = NULL;
2168: if (clPoints) *clPoints = NULL;
2169: }
2170: return(0);
2171: }
2175: /*@
2176: PetscSectionGetField - Get the subsection associated with a single field
2178: Input Parameters:
2179: + s - The PetscSection
2180: - field - The field number
2182: Output Parameter:
2183: . subs - The subsection for the given field
2185: Level: intermediate
2187: .seealso: PetscSectionSetNumFields()
2188: @*/
2189: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2190: {
2196: 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);
2197: PetscObjectReference((PetscObject) s->field[field]);
2198: *subs = s->field[field];
2199: return(0);
2200: }