Actual source code: vsectionis.c
petsc-3.6.4 2016-04-12
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: - point - the point
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
965: Output Parameter:
966: . gsection - The PetscSection for the global field layout
968: Note: This gives negative sizes and offsets to points not owned by this process
970: Level: developer
972: .seealso: PetscSectionCreate()
973: @*/
974: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscSection *gsection)
975: {
976: const PetscInt *pind = NULL;
977: PetscInt *recv = NULL, *neg = NULL;
978: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
979: PetscErrorCode ierr;
982: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
983: PetscSectionGetChart(s, &pStart, &pEnd);
984: PetscSectionSetChart(*gsection, pStart, pEnd);
985: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
986: nlocal = nroots; /* The local/leaf space matches global/root space */
987: /* Must allocate for all points visible to SF, which may be more than this section */
988: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
989: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
990: PetscMalloc2(nroots,&neg,nlocal,&recv);
991: PetscMemzero(neg,nroots*sizeof(PetscInt));
992: }
993: /* Mark all local points with negative dof */
994: for (p = pStart; p < pEnd; ++p) {
995: PetscSectionGetDof(s, p, &dof);
996: PetscSectionSetDof(*gsection, p, dof);
997: PetscSectionGetConstraintDof(s, p, &cdof);
998: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
999: if (neg) neg[p] = -(dof+1);
1000: }
1001: PetscSectionSetUpBC(*gsection);
1002: if (nroots >= 0) {
1003: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1004: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1005: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1006: for (p = pStart; p < pEnd; ++p) {
1007: if (recv[p] < 0) {
1008: (*gsection)->atlasDof[p-pStart] = recv[p];
1009: PetscSectionGetDof(s, p, &dof);
1010: 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);
1011: }
1012: }
1013: }
1014: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1015: if (s->perm) {ISGetIndices(s->perm, &pind);}
1016: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1017: const PetscInt q = pind ? pind[p] : p;
1019: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1020: (*gsection)->atlasOff[q] = off;
1021: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1022: }
1023: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1024: globalOff -= off;
1025: for (p = pStart, off = 0; p < pEnd; ++p) {
1026: (*gsection)->atlasOff[p-pStart] += globalOff;
1027: if (neg) neg[p] = -((*gsection)->atlasOff[p-pStart]+1);
1028: }
1029: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1030: /* Put in negative offsets for ghost points */
1031: if (nroots >= 0) {
1032: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1033: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1034: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1035: for (p = pStart; p < pEnd; ++p) {
1036: if (recv[p] < 0) (*gsection)->atlasOff[p-pStart] = recv[p];
1037: }
1038: }
1039: PetscFree2(neg,recv);
1040: PetscSectionViewFromOptions(*gsection,NULL,"-section_global_view");
1041: return(0);
1042: }
1046: /*@
1047: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1048: the local section and an SF describing the section point overlap.
1050: Input Parameters:
1051: + s - The PetscSection for the local field layout
1052: . sf - The SF describing parallel layout of the section points
1053: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1054: . numExcludes - The number of exclusion ranges
1055: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1057: Output Parameter:
1058: . gsection - The PetscSection for the global field layout
1060: Note: This gives negative sizes and offsets to points not owned by this process
1062: Level: developer
1064: .seealso: PetscSectionCreate()
1065: @*/
1066: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1067: {
1068: const PetscInt *pind = NULL;
1069: PetscInt *neg = NULL, *tmpOff = NULL;
1070: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1071: PetscErrorCode ierr;
1074: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1075: PetscSectionGetChart(s, &pStart, &pEnd);
1076: PetscSectionSetChart(*gsection, pStart, pEnd);
1077: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1078: if (nroots >= 0) {
1079: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1080: PetscCalloc1(nroots, &neg);
1081: if (nroots > pEnd-pStart) {
1082: PetscCalloc1(nroots, &tmpOff);
1083: } else {
1084: tmpOff = &(*gsection)->atlasDof[-pStart];
1085: }
1086: }
1087: /* Mark ghost points with negative dof */
1088: for (p = pStart; p < pEnd; ++p) {
1089: for (e = 0; e < numExcludes; ++e) {
1090: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1091: PetscSectionSetDof(*gsection, p, 0);
1092: break;
1093: }
1094: }
1095: if (e < numExcludes) continue;
1096: PetscSectionGetDof(s, p, &dof);
1097: PetscSectionSetDof(*gsection, p, dof);
1098: PetscSectionGetConstraintDof(s, p, &cdof);
1099: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1100: if (neg) neg[p] = -(dof+1);
1101: }
1102: PetscSectionSetUpBC(*gsection);
1103: if (nroots >= 0) {
1104: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1105: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1106: if (nroots > pEnd - pStart) {
1107: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1108: }
1109: }
1110: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1111: if (s->perm) {ISGetIndices(s->perm, &pind);}
1112: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1113: const PetscInt q = pind ? pind[p] : p;
1115: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1116: (*gsection)->atlasOff[q] = off;
1117: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1118: }
1119: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1120: globalOff -= off;
1121: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1122: (*gsection)->atlasOff[p] += globalOff;
1123: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1124: }
1125: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1126: /* Put in negative offsets for ghost points */
1127: if (nroots >= 0) {
1128: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1129: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1130: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1131: if (nroots > pEnd - pStart) {
1132: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1133: }
1134: }
1135: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1136: PetscFree(neg);
1137: return(0);
1138: }
1142: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1143: {
1144: PetscInt pStart, pEnd, p, localSize = 0;
1148: PetscSectionGetChart(s, &pStart, &pEnd);
1149: for (p = pStart; p < pEnd; ++p) {
1150: PetscInt dof;
1152: PetscSectionGetDof(s, p, &dof);
1153: if (dof > 0) ++localSize;
1154: }
1155: PetscLayoutCreate(comm, layout);
1156: PetscLayoutSetLocalSize(*layout, localSize);
1157: PetscLayoutSetBlockSize(*layout, 1);
1158: PetscLayoutSetUp(*layout);
1159: return(0);
1160: }
1164: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1165: {
1166: PetscInt pStart, pEnd, p, localSize = 0;
1170: PetscSectionGetChart(s, &pStart, &pEnd);
1171: for (p = pStart; p < pEnd; ++p) {
1172: PetscInt dof,cdof;
1174: PetscSectionGetDof(s, p, &dof);
1175: PetscSectionGetConstraintDof(s, p, &cdof);
1176: if (dof-cdof > 0) localSize += dof-cdof;
1177: }
1178: PetscLayoutCreate(comm, layout);
1179: PetscLayoutSetLocalSize(*layout, localSize);
1180: PetscLayoutSetBlockSize(*layout, 1);
1181: PetscLayoutSetUp(*layout);
1182: return(0);
1183: }
1187: /*@
1188: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1190: Not collective
1192: Input Parameters:
1193: + s - the PetscSection
1194: - point - the point
1196: Output Parameter:
1197: . offset - the offset
1199: Level: intermediate
1201: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1202: @*/
1203: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1204: {
1206: #if defined(PETSC_USE_DEBUG)
1207: 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);
1208: #endif
1209: *offset = s->atlasOff[point - s->pStart];
1210: return(0);
1211: }
1215: /*@
1216: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1218: Not collective
1220: Input Parameters:
1221: + s - the PetscSection
1222: . point - the point
1223: - offset - the offset
1225: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1227: Level: intermediate
1229: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1230: @*/
1231: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1232: {
1234: 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);
1235: s->atlasOff[point - s->pStart] = offset;
1236: return(0);
1237: }
1241: /*@
1242: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1244: Not collective
1246: Input Parameters:
1247: + s - the PetscSection
1248: . point - the point
1249: - field - the field
1251: Output Parameter:
1252: . offset - the offset
1254: Level: intermediate
1256: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1257: @*/
1258: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1259: {
1263: 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);
1264: PetscSectionGetOffset(s->field[field], point, offset);
1265: return(0);
1266: }
1270: /*@
1271: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1273: Not collective
1275: Input Parameters:
1276: + s - the PetscSection
1277: . point - the point
1278: . field - the field
1279: - offset - the offset
1281: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1283: Level: intermediate
1285: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1286: @*/
1287: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1288: {
1292: 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);
1293: PetscSectionSetOffset(s->field[field], point, offset);
1294: return(0);
1295: }
1299: /* This gives the offset on a point of the field, ignoring constraints */
1300: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1301: {
1302: PetscInt off, foff;
1306: 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);
1307: PetscSectionGetOffset(s, point, &off);
1308: PetscSectionGetOffset(s->field[field], point, &foff);
1309: *offset = foff - off;
1310: return(0);
1311: }
1315: /*@
1316: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1318: Not collective
1320: Input Parameter:
1321: . s - the PetscSection
1323: Output Parameters:
1324: + start - the minimum offset
1325: - end - one more than the maximum offset
1327: Level: intermediate
1329: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1330: @*/
1331: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1332: {
1333: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1337: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1338: PetscSectionGetChart(s, &pStart, &pEnd);
1339: for (p = 0; p < pEnd-pStart; ++p) {
1340: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1342: if (off >= 0) {
1343: os = PetscMin(os, off);
1344: oe = PetscMax(oe, off+dof);
1345: }
1346: }
1347: if (start) *start = os;
1348: if (end) *end = oe;
1349: return(0);
1350: }
1354: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, PetscInt fields[], PetscSection *subs)
1355: {
1356: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1360: if (!numFields) return(0);
1361: PetscSectionGetNumFields(s, &nF);
1362: if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1363: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1364: PetscSectionSetNumFields(*subs, numFields);
1365: for (f = 0; f < numFields; ++f) {
1366: const char *name = NULL;
1367: PetscInt numComp = 0;
1369: PetscSectionGetFieldName(s, fields[f], &name);
1370: PetscSectionSetFieldName(*subs, f, name);
1371: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1372: PetscSectionSetFieldComponents(*subs, f, numComp);
1373: }
1374: PetscSectionGetChart(s, &pStart, &pEnd);
1375: PetscSectionSetChart(*subs, pStart, pEnd);
1376: for (p = pStart; p < pEnd; ++p) {
1377: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1379: for (f = 0; f < numFields; ++f) {
1380: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1381: PetscSectionSetFieldDof(*subs, p, f, fdof);
1382: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1383: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1384: dof += fdof;
1385: cdof += cfdof;
1386: }
1387: PetscSectionSetDof(*subs, p, dof);
1388: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1389: maxCdof = PetscMax(cdof, maxCdof);
1390: }
1391: PetscSectionSetUp(*subs);
1392: if (maxCdof) {
1393: PetscInt *indices;
1395: PetscMalloc1(maxCdof, &indices);
1396: for (p = pStart; p < pEnd; ++p) {
1397: PetscInt cdof;
1399: PetscSectionGetConstraintDof(*subs, p, &cdof);
1400: if (cdof) {
1401: const PetscInt *oldIndices = NULL;
1402: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1404: for (f = 0; f < numFields; ++f) {
1405: PetscInt oldFoff = 0, oldf;
1407: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1408: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1409: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1410: /* This can be sped up if we assume sorted fields */
1411: for (oldf = 0; oldf < fields[f]; ++oldf) {
1412: PetscInt oldfdof = 0;
1413: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1414: oldFoff += oldfdof;
1415: }
1416: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1417: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1418: numConst += cfdof;
1419: fOff += fdof;
1420: }
1421: PetscSectionSetConstraintIndices(*subs, p, indices);
1422: }
1423: }
1424: PetscFree(indices);
1425: }
1426: return(0);
1427: }
1431: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1432: {
1433: const PetscInt *points = NULL, *indices = NULL;
1434: PetscInt numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;
1438: PetscSectionGetNumFields(s, &numFields);
1439: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1440: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1441: for (f = 0; f < numFields; ++f) {
1442: const char *name = NULL;
1443: PetscInt numComp = 0;
1445: PetscSectionGetFieldName(s, f, &name);
1446: PetscSectionSetFieldName(*subs, f, name);
1447: PetscSectionGetFieldComponents(s, f, &numComp);
1448: PetscSectionSetFieldComponents(*subs, f, numComp);
1449: }
1450: /* For right now, we do not try to squeeze the subchart */
1451: if (subpointMap) {
1452: ISGetSize(subpointMap, &numSubpoints);
1453: ISGetIndices(subpointMap, &points);
1454: }
1455: PetscSectionGetChart(s, &pStart, &pEnd);
1456: PetscSectionSetChart(*subs, 0, numSubpoints);
1457: for (p = pStart; p < pEnd; ++p) {
1458: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1460: PetscFindInt(p, numSubpoints, points, &subp);
1461: if (subp < 0) continue;
1462: for (f = 0; f < numFields; ++f) {
1463: PetscSectionGetFieldDof(s, p, f, &fdof);
1464: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1465: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1466: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1467: }
1468: PetscSectionGetDof(s, p, &dof);
1469: PetscSectionSetDof(*subs, subp, dof);
1470: PetscSectionGetConstraintDof(s, p, &cdof);
1471: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1472: }
1473: PetscSectionSetUp(*subs);
1474: /* Change offsets to original offsets */
1475: for (p = pStart; p < pEnd; ++p) {
1476: PetscInt off, foff = 0;
1478: PetscFindInt(p, numSubpoints, points, &subp);
1479: if (subp < 0) continue;
1480: for (f = 0; f < numFields; ++f) {
1481: PetscSectionGetFieldOffset(s, p, f, &foff);
1482: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1483: }
1484: PetscSectionGetOffset(s, p, &off);
1485: PetscSectionSetOffset(*subs, subp, off);
1486: }
1487: /* Copy constraint indices */
1488: for (subp = 0; subp < numSubpoints; ++subp) {
1489: PetscInt cdof;
1491: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1492: if (cdof) {
1493: for (f = 0; f < numFields; ++f) {
1494: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1495: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1496: }
1497: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1498: PetscSectionSetConstraintIndices(*subs, subp, indices);
1499: }
1500: }
1501: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1502: return(0);
1503: }
1507: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1508: {
1509: PetscInt p;
1510: PetscMPIInt rank;
1514: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1515: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
1516: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1517: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1518: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1519: PetscInt b;
1521: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1522: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1523: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1524: }
1525: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1526: } else {
1527: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1528: }
1529: }
1530: PetscViewerFlush(viewer);
1531: return(0);
1532: }
1536: /*@C
1537: PetscSectionView - Views a PetscSection
1539: Collective on PetscSection
1541: Input Parameters:
1542: + s - the PetscSection object to view
1543: - v - the viewer
1545: Level: developer
1547: .seealso PetscSectionCreate(), PetscSectionDestroy()
1548: @*/
1549: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1550: {
1551: PetscBool isascii;
1552: PetscInt f;
1556: if (!viewer) {PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer);}
1558: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1559: if (isascii) {
1560: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1561: if (s->numFields) {
1562: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1563: for (f = 0; f < s->numFields; ++f) {
1564: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
1565: PetscSectionView_ASCII(s->field[f], viewer);
1566: }
1567: } else {
1568: PetscSectionView_ASCII(s, viewer);
1569: }
1570: }
1571: return(0);
1572: }
1576: /*@
1577: PetscSectionReset - Frees all section data.
1579: Not collective
1581: Input Parameters:
1582: . s - the PetscSection
1584: Level: developer
1586: .seealso: PetscSection, PetscSectionCreate()
1587: @*/
1588: PetscErrorCode PetscSectionReset(PetscSection s)
1589: {
1590: PetscInt f;
1594: PetscFree(s->numFieldComponents);
1595: for (f = 0; f < s->numFields; ++f) {
1596: PetscSectionDestroy(&s->field[f]);
1597: PetscFree(s->fieldNames[f]);
1598: }
1599: PetscFree(s->fieldNames);
1600: PetscFree(s->field);
1601: PetscSectionDestroy(&s->bc);
1602: PetscFree(s->bcIndices);
1603: PetscFree2(s->atlasDof, s->atlasOff);
1604: PetscSectionDestroy(&s->clSection);
1605: ISDestroy(&s->clPoints);
1606: ISDestroy(&s->perm);
1608: s->pStart = -1;
1609: s->pEnd = -1;
1610: s->maxDof = 0;
1611: s->setup = PETSC_FALSE;
1612: s->numFields = 0;
1613: s->clObj = NULL;
1614: return(0);
1615: }
1619: /*@
1620: PetscSectionDestroy - Frees a section object and frees its range if that exists.
1622: Not collective
1624: Input Parameters:
1625: . s - the PetscSection
1627: Level: developer
1629: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1630: recommended they not be used in user codes unless you really gain something in their use.
1632: .seealso: PetscSection, PetscSectionCreate()
1633: @*/
1634: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1635: {
1639: if (!*s) return(0);
1641: if (--((PetscObject)(*s))->refct > 0) {
1642: *s = NULL;
1643: return(0);
1644: }
1645: PetscSectionReset(*s);
1646: PetscHeaderDestroy(s);
1647: return(0);
1648: }
1652: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1653: {
1654: const PetscInt p = point - s->pStart;
1657: *values = &baseArray[s->atlasOff[p]];
1658: return(0);
1659: }
1663: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1664: {
1665: PetscInt *array;
1666: const PetscInt p = point - s->pStart;
1667: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1668: PetscInt cDim = 0;
1672: PetscSectionGetConstraintDof(s, p, &cDim);
1673: array = &baseArray[s->atlasOff[p]];
1674: if (!cDim) {
1675: if (orientation >= 0) {
1676: const PetscInt dim = s->atlasDof[p];
1677: PetscInt i;
1679: if (mode == INSERT_VALUES) {
1680: for (i = 0; i < dim; ++i) array[i] = values[i];
1681: } else {
1682: for (i = 0; i < dim; ++i) array[i] += values[i];
1683: }
1684: } else {
1685: PetscInt offset = 0;
1686: PetscInt j = -1, field, i;
1688: for (field = 0; field < s->numFields; ++field) {
1689: const PetscInt dim = s->field[field]->atlasDof[p];
1691: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1692: offset += dim;
1693: }
1694: }
1695: } else {
1696: if (orientation >= 0) {
1697: const PetscInt dim = s->atlasDof[p];
1698: PetscInt cInd = 0, i;
1699: const PetscInt *cDof;
1701: PetscSectionGetConstraintIndices(s, point, &cDof);
1702: if (mode == INSERT_VALUES) {
1703: for (i = 0; i < dim; ++i) {
1704: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1705: array[i] = values[i];
1706: }
1707: } else {
1708: for (i = 0; i < dim; ++i) {
1709: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1710: array[i] += values[i];
1711: }
1712: }
1713: } else {
1714: const PetscInt *cDof;
1715: PetscInt offset = 0;
1716: PetscInt cOffset = 0;
1717: PetscInt j = 0, field;
1719: PetscSectionGetConstraintIndices(s, point, &cDof);
1720: for (field = 0; field < s->numFields; ++field) {
1721: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
1722: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1723: const PetscInt sDim = dim - tDim;
1724: PetscInt cInd = 0, i,k;
1726: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1727: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1728: array[j] = values[k];
1729: }
1730: offset += dim;
1731: cOffset += dim - tDim;
1732: }
1733: }
1734: }
1735: return(0);
1736: }
1740: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1741: {
1745: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1746: return(0);
1747: }
1751: /*@C
1752: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
1754: Input Parameters:
1755: + s - The PetscSection
1756: - point - The point
1758: Output Parameter:
1759: . indices - The constrained dofs
1761: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
1763: Level: advanced
1765: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1766: @*/
1767: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1768: {
1772: if (s->bc) {
1773: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1774: } else *indices = NULL;
1775: return(0);
1776: }
1780: /*@C
1781: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
1783: Input Parameters:
1784: + s - The PetscSection
1785: . point - The point
1786: - indices - The constrained dofs
1788: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
1790: Level: advanced
1792: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1793: @*/
1794: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1795: {
1799: if (s->bc) {
1800: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1801: }
1802: return(0);
1803: }
1807: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1808: {
1812: 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);
1813: PetscSectionGetConstraintIndices(s->field[field], point, indices);
1814: return(0);
1815: }
1819: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
1820: {
1824: 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);
1825: PetscSectionSetConstraintIndices(s->field[field], point, indices);
1826: return(0);
1827: }
1831: /*@
1832: PetscSectionPermute - Reorder the section according to the input point permutation
1834: Collective on PetscSection
1836: Input Parameter:
1837: + section - The PetscSection object
1838: - perm - The point permutation, old point p becomes new point perm[p]
1840: Output Parameter:
1841: . sectionNew - The permuted PetscSection
1843: Level: intermediate
1845: .keywords: mesh
1846: .seealso: MatPermute()
1847: @*/
1848: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1849: {
1850: PetscSection s = section, sNew;
1851: const PetscInt *perm;
1852: PetscInt numFields, f, numPoints, pStart, pEnd, p;
1853: PetscErrorCode ierr;
1859: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1860: PetscSectionGetNumFields(s, &numFields);
1861: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1862: for (f = 0; f < numFields; ++f) {
1863: const char *name;
1864: PetscInt numComp;
1866: PetscSectionGetFieldName(s, f, &name);
1867: PetscSectionSetFieldName(sNew, f, name);
1868: PetscSectionGetFieldComponents(s, f, &numComp);
1869: PetscSectionSetFieldComponents(sNew, f, numComp);
1870: }
1871: ISGetLocalSize(permutation, &numPoints);
1872: ISGetIndices(permutation, &perm);
1873: PetscSectionGetChart(s, &pStart, &pEnd);
1874: PetscSectionSetChart(sNew, pStart, pEnd);
1875: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1876: for (p = pStart; p < pEnd; ++p) {
1877: PetscInt dof, cdof;
1879: PetscSectionGetDof(s, p, &dof);
1880: PetscSectionSetDof(sNew, perm[p], dof);
1881: PetscSectionGetConstraintDof(s, p, &cdof);
1882: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1883: for (f = 0; f < numFields; ++f) {
1884: PetscSectionGetFieldDof(s, p, f, &dof);
1885: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1886: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1887: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1888: }
1889: }
1890: PetscSectionSetUp(sNew);
1891: for (p = pStart; p < pEnd; ++p) {
1892: const PetscInt *cind;
1893: PetscInt cdof;
1895: PetscSectionGetConstraintDof(s, p, &cdof);
1896: if (cdof) {
1897: PetscSectionGetConstraintIndices(s, p, &cind);
1898: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1899: }
1900: for (f = 0; f < numFields; ++f) {
1901: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1902: if (cdof) {
1903: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1904: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1905: }
1906: }
1907: }
1908: ISRestoreIndices(permutation, &perm);
1909: *sectionNew = sNew;
1910: return(0);
1911: }
1915: /*@C
1916: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
1918: Collective
1920: Input Parameters:
1921: + sf - The SF
1922: - rootSection - Section defined on root space
1924: Output Parameters:
1925: + remoteOffsets - root offsets in leaf storage, or NULL
1926: - leafSection - Section defined on the leaf space
1928: Level: intermediate
1930: .seealso: PetscSFCreate()
1931: @*/
1932: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1933: {
1934: PetscSF embedSF;
1935: const PetscInt *ilocal, *indices;
1936: IS selected;
1937: PetscInt numFields, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;
1941: PetscSectionGetNumFields(rootSection, &numFields);
1942: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
1943: for (f = 0; f < numFields; ++f) {
1944: PetscInt numComp = 0;
1945: PetscSectionGetFieldComponents(rootSection, f, &numComp);
1946: PetscSectionSetFieldComponents(leafSection, f, numComp);
1947: }
1948: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1949: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1950: ISGetIndices(selected, &indices);
1951: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1952: ISRestoreIndices(selected, &indices);
1953: ISDestroy(&selected);
1954: PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
1955: if (nleaves && ilocal) {
1956: for (i = 0; i < nleaves; ++i) {
1957: lpStart = PetscMin(lpStart, ilocal[i]);
1958: lpEnd = PetscMax(lpEnd, ilocal[i]);
1959: }
1960: ++lpEnd;
1961: } else {
1962: lpStart = 0;
1963: lpEnd = nleaves;
1964: }
1965: PetscSectionSetChart(leafSection, lpStart, lpEnd);
1966: /* Could fuse these at the cost of a copy and extra allocation */
1967: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1968: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1969: for (f = 0; f < numFields; ++f) {
1970: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1971: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1972: }
1973: if (remoteOffsets) {
1974: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
1975: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1976: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1977: }
1978: PetscSFDestroy(&embedSF);
1979: PetscSectionSetUp(leafSection);
1980: return(0);
1981: }
1985: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
1986: {
1987: PetscSF embedSF;
1988: const PetscInt *indices;
1989: IS selected;
1990: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
1991: PetscErrorCode ierr;
1994: *remoteOffsets = NULL;
1995: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
1996: if (numRoots < 0) return(0);
1997: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1998: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1999: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2000: ISGetIndices(selected, &indices);
2001: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2002: ISRestoreIndices(selected, &indices);
2003: ISDestroy(&selected);
2004: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2005: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2006: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2007: PetscSFDestroy(&embedSF);
2008: return(0);
2009: }
2013: /*@C
2014: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2016: Input Parameters:
2017: + sf - The SF
2018: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section), or NULL
2019: - remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2021: Output Parameters:
2022: + leafSection - Data layout of local points for incoming data (this is the distributed section)
2023: - sectionSF - The new SF
2025: Note: Either rootSection or remoteOffsets can be specified
2027: Level: intermediate
2029: .seealso: PetscSFCreate()
2030: @*/
2031: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2032: {
2033: MPI_Comm comm;
2034: const PetscInt *localPoints;
2035: const PetscSFNode *remotePoints;
2036: PetscInt lpStart, lpEnd;
2037: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2038: PetscInt *localIndices;
2039: PetscSFNode *remoteIndices;
2040: PetscInt i, ind;
2041: PetscErrorCode ierr;
2049: PetscObjectGetComm((PetscObject)sf,&comm);
2050: PetscSFCreate(comm, sectionSF);
2051: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2052: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2053: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2054: if (numRoots < 0) return(0);
2055: for (i = 0; i < numPoints; ++i) {
2056: PetscInt localPoint = localPoints ? localPoints[i] : i;
2057: PetscInt dof;
2059: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2060: PetscSectionGetDof(leafSection, localPoint, &dof);
2061: numIndices += dof;
2062: }
2063: }
2064: PetscMalloc1(numIndices, &localIndices);
2065: PetscMalloc1(numIndices, &remoteIndices);
2066: /* Create new index graph */
2067: for (i = 0, ind = 0; i < numPoints; ++i) {
2068: PetscInt localPoint = localPoints ? localPoints[i] : i;
2069: PetscInt rank = remotePoints[i].rank;
2071: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2072: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2073: PetscInt loff, dof, d;
2075: PetscSectionGetOffset(leafSection, localPoint, &loff);
2076: PetscSectionGetDof(leafSection, localPoint, &dof);
2077: for (d = 0; d < dof; ++d, ++ind) {
2078: localIndices[ind] = loff+d;
2079: remoteIndices[ind].rank = rank;
2080: remoteIndices[ind].index = remoteOffset+d;
2081: }
2082: }
2083: }
2084: PetscFree(remoteOffsets);
2085: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2086: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2087: return(0);
2088: }
2092: /*@
2093: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2095: Input Parameters:
2096: + section - The PetscSection
2097: . obj - A PetscObject which serves as the key for this index
2098: . clSection - Section giving the size of the closure of each point
2099: - clPoints - IS giving the points in each closure
2101: Note: We compress out closure points with no dofs in this section
2103: Level: intermediate
2105: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2106: @*/
2107: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2108: {
2112: section->clObj = obj;
2113: PetscSectionDestroy(§ion->clSection);
2114: ISDestroy(§ion->clPoints);
2115: section->clSection = clSection;
2116: section->clPoints = clPoints;
2117: return(0);
2118: }
2122: /*@
2123: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2125: Input Parameters:
2126: + section - The PetscSection
2127: - obj - A PetscObject which serves as the key for this index
2129: Output Parameters:
2130: + clSection - Section giving the size of the closure of each point
2131: - clPoints - IS giving the points in each closure
2133: Note: We compress out closure points with no dofs in this section
2135: Level: intermediate
2137: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2138: @*/
2139: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2140: {
2142: if (section->clObj == obj) {
2143: if (clSection) *clSection = section->clSection;
2144: if (clPoints) *clPoints = section->clPoints;
2145: } else {
2146: if (clSection) *clSection = NULL;
2147: if (clPoints) *clPoints = NULL;
2148: }
2149: return(0);
2150: }
2154: /*@
2155: PetscSectionGetField - Get the subsection associated with a single field
2157: Input Parameters:
2158: + s - The PetscSection
2159: - field - The field number
2161: Output Parameter:
2162: . subs - The subsection for the given field
2164: Level: intermediate
2166: .seealso: PetscSectionSetNumFields()
2167: @*/
2168: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2169: {
2175: 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);
2176: PetscObjectReference((PetscObject) s->field[field]);
2177: *subs = s->field[field];
2178: return(0);
2179: }