Actual source code: vsectionis.c
petsc-3.5.4 2015-05-23
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,_p_PetscSection,int,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: PetscSectionClone - 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 PetscSectionClone(PetscSection section, PetscSection *newSection)
84: {
85: IS perm;
86: PetscInt numFields, f, pStart, pEnd, p;
90: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
91: PetscSectionGetNumFields(section, &numFields);
92: if (numFields) {PetscSectionSetNumFields(*newSection, numFields);}
93: for (f = 0; f < numFields; ++f) {
94: const char *name = NULL;
95: PetscInt numComp = 0;
97: PetscSectionGetFieldName(section, f, &name);
98: PetscSectionSetFieldName(*newSection, f, name);
99: PetscSectionGetFieldComponents(section, f, &numComp);
100: PetscSectionSetFieldComponents(*newSection, f, numComp);
101: }
102: PetscSectionGetChart(section, &pStart, &pEnd);
103: PetscSectionSetChart(*newSection, pStart, pEnd);
104: PetscSectionGetPermutation(section, &perm);
105: PetscSectionSetPermutation(*newSection, perm);
106: for (p = pStart; p < pEnd; ++p) {
107: PetscInt dof, cdof, fcdof = 0;
109: PetscSectionGetDof(section, p, &dof);
110: PetscSectionSetDof(*newSection, p, dof);
111: PetscSectionGetConstraintDof(section, p, &cdof);
112: if (cdof) {PetscSectionSetConstraintDof(*newSection, p, cdof);}
113: for (f = 0; f < numFields; ++f) {
114: PetscSectionGetFieldDof(section, p, f, &dof);
115: PetscSectionSetFieldDof(*newSection, p, f, dof);
116: if (cdof) {
117: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
118: if (fcdof) {PetscSectionSetFieldConstraintDof(*newSection, p, f, fcdof);}
119: }
120: }
121: }
122: PetscSectionSetUp(*newSection);
123: for (p = pStart; p < pEnd; ++p) {
124: PetscInt off, cdof, fcdof = 0;
125: const PetscInt *cInd;
127: /* Must set offsets in case they do not agree with the prefix sums */
128: PetscSectionGetOffset(section, p, &off);
129: PetscSectionSetOffset(*newSection, p, off);
130: PetscSectionGetConstraintDof(section, p, &cdof);
131: if (cdof) {
132: PetscSectionGetConstraintIndices(section, p, &cInd);
133: PetscSectionSetConstraintIndices(*newSection, p, cInd);
134: for (f = 0; f < numFields; ++f) {
135: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
136: if (fcdof) {
137: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
138: PetscSectionSetFieldConstraintIndices(*newSection, p, f, cInd);
139: }
140: }
141: }
142: }
143: return(0);
144: }
148: /*@
149: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
151: Not collective
153: Input Parameter:
154: . s - the PetscSection
156: Output Parameter:
157: . numFields - the number of fields defined, or 0 if none were defined
159: Level: intermediate
161: .seealso: PetscSectionSetNumFields()
162: @*/
163: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
164: {
167: *numFields = s->numFields;
168: return(0);
169: }
173: /*@
174: PetscSectionSetNumFields - Sets the number of fields.
176: Not collective
178: Input Parameters:
179: + s - the PetscSection
180: - numFields - the number of fields
182: Level: intermediate
184: .seealso: PetscSectionGetNumFields()
185: @*/
186: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
187: {
188: PetscInt f;
192: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %d must be positive", numFields);
193: PetscSectionReset(s);
195: s->numFields = numFields;
196: PetscMalloc1(s->numFields, &s->numFieldComponents);
197: PetscMalloc1(s->numFields, &s->fieldNames);
198: PetscMalloc1(s->numFields, &s->field);
199: for (f = 0; f < s->numFields; ++f) {
200: char name[64];
202: s->numFieldComponents[f] = 1;
204: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
205: PetscSNPrintf(name, 64, "Field_%D", f);
206: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
207: }
208: return(0);
209: }
213: /*@C
214: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
216: Not Collective
218: Input Parameters:
219: + s - the PetscSection
220: - field - the field number
222: Output Parameter:
223: . fieldName - the field name
225: Level: developer
227: .seealso: PetscSectionSetFieldName()
228: @*/
229: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
230: {
233: 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);
234: *fieldName = s->fieldNames[field];
235: return(0);
236: }
240: /*@C
241: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
243: Not Collective
245: Input Parameters:
246: + s - the PetscSection
247: . field - the field number
248: - fieldName - the field name
250: Level: developer
252: .seealso: PetscSectionGetFieldName()
253: @*/
254: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
255: {
260: 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);
261: PetscFree(s->fieldNames[field]);
262: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
263: return(0);
264: }
268: /*@
269: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
271: Not collective
273: Input Parameters:
274: + s - the PetscSection
275: - field - the field number
277: Output Parameter:
278: . numComp - the number of field components
280: Level: intermediate
282: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
283: @*/
284: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
285: {
288: 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);
289: *numComp = s->numFieldComponents[field];
290: return(0);
291: }
295: /*@
296: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
298: Not collective
300: Input Parameters:
301: + s - the PetscSection
302: . field - the field number
303: - numComp - the number of field components
305: Level: intermediate
307: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
308: @*/
309: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
310: {
312: 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);
313: s->numFieldComponents[field] = numComp;
314: return(0);
315: }
319: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
320: {
324: if (!s->bc) {
325: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
326: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
327: }
328: return(0);
329: }
333: /*@
334: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points in the lie.
336: Not collective
338: Input Parameter:
339: . s - the PetscSection
341: Output Parameters:
342: + pStart - the first point
343: - pEnd - one past the last point
345: Level: intermediate
347: .seealso: PetscSectionSetChart(), PetscSectionCreate()
348: @*/
349: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
350: {
353: if (pStart) *pStart = s->pStart;
354: if (pEnd) *pEnd = s->pEnd;
355: return(0);
356: }
360: /*@
361: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points in the lie.
363: Not collective
365: Input Parameters:
366: + s - the PetscSection
367: . pStart - the first point
368: - pEnd - one past the last point
370: Level: intermediate
372: .seealso: PetscSectionGetChart(), PetscSectionCreate()
373: @*/
374: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
375: {
376: PetscInt f;
381: /* Cannot Reset() because it destroys field information */
382: s->setup = PETSC_FALSE;
383: PetscSectionDestroy(&s->bc);
384: PetscFree(s->bcIndices);
385: PetscFree2(s->atlasDof, s->atlasOff);
387: s->pStart = pStart;
388: s->pEnd = pEnd;
389: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
390: PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
391: for (f = 0; f < s->numFields; ++f) {
392: PetscSectionSetChart(s->field[f], pStart, pEnd);
393: }
394: return(0);
395: }
399: /*@
400: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
402: Not collective
404: Input Parameter:
405: . s - the PetscSection
407: Output Parameters:
408: . perm - The permutation as an IS
410: Level: intermediate
412: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
413: @*/
414: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
415: {
419: return(0);
420: }
424: /*@
425: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
427: Not collective
429: Input Parameters:
430: + s - the PetscSection
431: - perm - the permutation of points
433: Level: intermediate
435: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
436: @*/
437: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
438: {
442: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
443: if (s->perm != perm) {
444: ISDestroy(&s->perm);
445: s->perm = perm;
446: PetscObjectReference((PetscObject) s->perm);
447: }
448: return(0);
449: }
453: /*@
454: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
456: Not collective
458: Input Parameters:
459: + s - the PetscSection
460: - point - the point
462: Output Parameter:
463: . numDof - the number of dof
465: Level: intermediate
467: .seealso: PetscSectionSetDof(), PetscSectionCreate()
468: @*/
469: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
470: {
472: #if defined(PETSC_USE_DEBUG)
473: 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);
474: #endif
475: *numDof = s->atlasDof[point - s->pStart];
476: return(0);
477: }
481: /*@
482: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
484: Not collective
486: Input Parameters:
487: + s - the PetscSection
488: . point - the point
489: - numDof - the number of dof
491: Level: intermediate
493: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
494: @*/
495: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
496: {
498: 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);
499: s->atlasDof[point - s->pStart] = numDof;
500: return(0);
501: }
505: /*@
506: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
508: Not collective
510: Input Parameters:
511: + s - the PetscSection
512: . point - the point
513: - numDof - the number of additional dof
515: Level: intermediate
517: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
518: @*/
519: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
520: {
522: 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);
523: s->atlasDof[point - s->pStart] += numDof;
524: return(0);
525: }
529: /*@
530: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
532: Not collective
534: Input Parameters:
535: + s - the PetscSection
536: . point - the point
537: - field - the field
539: Output Parameter:
540: . numDof - the number of dof
542: Level: intermediate
544: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
545: @*/
546: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
547: {
551: 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);
552: PetscSectionGetDof(s->field[field], point, numDof);
553: return(0);
554: }
558: /*@
559: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
561: Not collective
563: Input Parameters:
564: + s - the PetscSection
565: . point - the point
566: . field - the field
567: - numDof - the number of dof
569: Level: intermediate
571: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
572: @*/
573: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
574: {
578: 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);
579: PetscSectionSetDof(s->field[field], point, numDof);
580: return(0);
581: }
585: /*@
586: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
588: Not collective
590: Input Parameters:
591: + s - the PetscSection
592: . point - the point
593: . field - the field
594: - numDof - the number of dof
596: Level: intermediate
598: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
599: @*/
600: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
601: {
605: 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);
606: PetscSectionAddDof(s->field[field], point, numDof);
607: return(0);
608: }
612: /*@
613: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
615: Not collective
617: Input Parameters:
618: + s - the PetscSection
619: - point - the point
621: Output Parameter:
622: . numDof - the number of dof which are fixed by constraints
624: Level: intermediate
626: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
627: @*/
628: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
629: {
633: if (s->bc) {
634: PetscSectionGetDof(s->bc, point, numDof);
635: } else *numDof = 0;
636: return(0);
637: }
641: /*@
642: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
644: Not collective
646: Input Parameters:
647: + s - the PetscSection
648: . point - the point
649: - numDof - the number of dof which are fixed by constraints
651: Level: intermediate
653: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
654: @*/
655: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
656: {
660: if (numDof) {
661: PetscSectionCheckConstraints_Static(s);
662: PetscSectionSetDof(s->bc, point, numDof);
663: }
664: return(0);
665: }
669: /*@
670: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
672: Not collective
674: Input Parameters:
675: + s - the PetscSection
676: . point - the point
677: - numDof - the number of additional dof which are fixed by constraints
679: Level: intermediate
681: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
682: @*/
683: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
684: {
688: if (numDof) {
689: PetscSectionCheckConstraints_Static(s);
690: PetscSectionAddDof(s->bc, point, numDof);
691: }
692: return(0);
693: }
697: /*@
698: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
700: Not collective
702: Input Parameters:
703: + s - the PetscSection
704: . point - the point
705: - field - the field
707: Output Parameter:
708: . numDof - the number of dof which are fixed by constraints
710: Level: intermediate
712: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
713: @*/
714: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
715: {
719: 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);
720: PetscSectionGetConstraintDof(s->field[field], point, numDof);
721: return(0);
722: }
726: /*@
727: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
729: Not collective
731: Input Parameters:
732: + s - the PetscSection
733: . point - the point
734: . field - the field
735: - numDof - the number of dof which are fixed by constraints
737: Level: intermediate
739: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
740: @*/
741: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
742: {
746: 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);
747: PetscSectionSetConstraintDof(s->field[field], point, numDof);
748: return(0);
749: }
753: /*@
754: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
756: Not collective
758: Input Parameters:
759: + s - the PetscSection
760: . point - the point
761: . field - the field
762: - numDof - the number of additional dof which are fixed by constraints
764: Level: intermediate
766: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
767: @*/
768: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
769: {
773: 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);
774: PetscSectionAddConstraintDof(s->field[field], point, numDof);
775: return(0);
776: }
780: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
781: {
785: if (s->bc) {
786: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
788: PetscSectionSetUp(s->bc);
789: PetscMalloc1((s->bc->atlasOff[last] + s->bc->atlasDof[last]), &s->bcIndices);
790: }
791: return(0);
792: }
796: /*@
797: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
799: Not collective
801: Input Parameter:
802: . s - the PetscSection
804: Level: intermediate
806: .seealso: PetscSectionCreate()
807: @*/
808: PetscErrorCode PetscSectionSetUp(PetscSection s)
809: {
810: const PetscInt *pind = NULL;
811: PetscInt offset = 0, p, f;
812: PetscErrorCode ierr;
815: if (s->setup) return(0);
816: s->setup = PETSC_TRUE;
817: if (s->perm) {ISGetIndices(s->perm, &pind);}
818: for (p = 0; p < s->pEnd - s->pStart; ++p) {
819: const PetscInt q = pind ? pind[p] : p;
821: s->atlasOff[q] = offset;
822: offset += s->atlasDof[q];
823: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
824: }
825: PetscSectionSetUpBC(s);
826: /* Assume that all fields have the same chart */
827: for (p = 0; p < s->pEnd - s->pStart; ++p) {
828: const PetscInt q = pind ? pind[p] : p;
829: PetscInt off = s->atlasOff[q];
831: for (f = 0; f < s->numFields; ++f) {
832: PetscSection sf = s->field[f];
834: sf->atlasOff[q] = off;
835: off += sf->atlasDof[q];
836: }
837: }
838: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
839: for (f = 0; f < s->numFields; ++f) {
840: PetscSectionSetUpBC(s->field[f]);
841: }
842: return(0);
843: }
847: /*@
848: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
850: Not collective
852: Input Parameters:
853: . s - the PetscSection
855: Output Parameter:
856: . maxDof - the maximum dof
858: Level: intermediate
860: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
861: @*/
862: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
863: {
865: *maxDof = s->maxDof;
866: return(0);
867: }
871: /*@
872: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
874: Not collective
876: Input Parameters:
877: + s - the PetscSection
878: - point - the point
880: Output Parameter:
881: . size - the size of an array which can hold all the dofs
883: Level: intermediate
885: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
886: @*/
887: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
888: {
889: PetscInt p, n = 0;
892: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
893: *size = n;
894: return(0);
895: }
899: /*@
900: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
902: Not collective
904: Input Parameters:
905: + s - the PetscSection
906: - point - the point
908: Output Parameter:
909: . size - the size of an array which can hold all unconstrained dofs
911: Level: intermediate
913: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
914: @*/
915: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
916: {
917: PetscInt p, n = 0;
920: for (p = 0; p < s->pEnd - s->pStart; ++p) {
921: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
922: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
923: }
924: *size = n;
925: return(0);
926: }
930: /*@
931: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
932: the local section and an SF describing the section point overlap.
934: Input Parameters:
935: + s - The PetscSection for the local field layout
936: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
937: - includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
939: Output Parameter:
940: . gsection - The PetscSection for the global field layout
942: Note: This gives negative sizes and offsets to points not owned by this process
944: Level: developer
946: .seealso: PetscSectionCreate()
947: @*/
948: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscSection *gsection)
949: {
950: const PetscInt *pind = NULL;
951: PetscInt *recv = NULL, *neg = NULL;
952: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
953: PetscErrorCode ierr;
956: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
957: PetscSectionGetChart(s, &pStart, &pEnd);
958: PetscSectionSetChart(*gsection, pStart, pEnd);
959: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
960: nlocal = nroots; /* The local/leaf space matches global/root space */
961: /* Must allocate for all points visible to SF, which may be more than this section */
962: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
963: PetscMalloc2(nroots,&neg,nlocal,&recv);
964: PetscMemzero(neg,nroots*sizeof(PetscInt));
965: }
966: /* Mark all local points with negative dof */
967: for (p = pStart; p < pEnd; ++p) {
968: PetscSectionGetDof(s, p, &dof);
969: PetscSectionSetDof(*gsection, p, dof);
970: PetscSectionGetConstraintDof(s, p, &cdof);
971: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
972: if (neg) neg[p] = -(dof+1);
973: }
974: PetscSectionSetUpBC(*gsection);
975: if (nroots >= 0) {
976: PetscMemzero(recv,nlocal*sizeof(PetscInt));
977: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
978: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
979: for (p = pStart; p < pEnd; ++p) {
980: if (recv[p] < 0) {
981: (*gsection)->atlasDof[p-pStart] = recv[p];
982: PetscSectionGetDof(s, p, &dof);
983: 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);
984: }
985: }
986: }
987: /* Calculate new sizes, get proccess offset, and calculate point offsets */
988: if (s->perm) {ISGetIndices(s->perm, &pind);}
989: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
990: const PetscInt q = pind ? pind[p] : p;
992: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
993: (*gsection)->atlasOff[q] = off;
994: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
995: }
996: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
997: globalOff -= off;
998: for (p = pStart, off = 0; p < pEnd; ++p) {
999: (*gsection)->atlasOff[p-pStart] += globalOff;
1000: if (neg) neg[p] = -((*gsection)->atlasOff[p-pStart]+1);
1001: }
1002: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1003: /* Put in negative offsets for ghost points */
1004: if (nroots >= 0) {
1005: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1006: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1007: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1008: for (p = pStart; p < pEnd; ++p) {
1009: if (recv[p] < 0) (*gsection)->atlasOff[p-pStart] = recv[p];
1010: }
1011: }
1012: PetscFree2(neg,recv);
1013: PetscSectionViewFromOptions(*gsection,NULL,"-section_global_view");
1014: return(0);
1015: }
1019: /*@
1020: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1021: the local section and an SF describing the section point overlap.
1023: Input Parameters:
1024: + s - The PetscSection for the local field layout
1025: . sf - The SF describing parallel layout of the section points
1026: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1027: . numExcludes - The number of exclusion ranges
1028: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1030: Output Parameter:
1031: . gsection - The PetscSection for the global field layout
1033: Note: This gives negative sizes and offsets to points not owned by this process
1035: Level: developer
1037: .seealso: PetscSectionCreate()
1038: @*/
1039: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1040: {
1041: const PetscInt *pind = NULL;
1042: PetscInt *neg = NULL, *tmpOff = NULL;
1043: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1044: PetscErrorCode ierr;
1047: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1048: PetscSectionGetChart(s, &pStart, &pEnd);
1049: PetscSectionSetChart(*gsection, pStart, pEnd);
1050: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1051: if (nroots >= 0) {
1052: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1053: PetscCalloc1(nroots, &neg);
1054: if (nroots > pEnd-pStart) {
1055: PetscCalloc1(nroots, &tmpOff);
1056: } else {
1057: tmpOff = &(*gsection)->atlasDof[-pStart];
1058: }
1059: }
1060: /* Mark ghost points with negative dof */
1061: for (p = pStart; p < pEnd; ++p) {
1062: for (e = 0; e < numExcludes; ++e) {
1063: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1064: PetscSectionSetDof(*gsection, p, 0);
1065: break;
1066: }
1067: }
1068: if (e < numExcludes) continue;
1069: PetscSectionGetDof(s, p, &dof);
1070: PetscSectionSetDof(*gsection, p, dof);
1071: PetscSectionGetConstraintDof(s, p, &cdof);
1072: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1073: if (neg) neg[p] = -(dof+1);
1074: }
1075: PetscSectionSetUpBC(*gsection);
1076: if (nroots >= 0) {
1077: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1078: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1079: if (nroots > pEnd - pStart) {
1080: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1081: }
1082: }
1083: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1084: if (s->perm) {ISGetIndices(s->perm, &pind);}
1085: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1086: const PetscInt q = pind ? pind[p] : p;
1088: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1089: (*gsection)->atlasOff[q] = off;
1090: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1091: }
1092: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1093: globalOff -= off;
1094: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1095: (*gsection)->atlasOff[p] += globalOff;
1096: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1097: }
1098: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1099: /* Put in negative offsets for ghost points */
1100: if (nroots >= 0) {
1101: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1102: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1103: if (nroots > pEnd - pStart) {
1104: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1105: }
1106: }
1107: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1108: PetscFree(neg);
1109: return(0);
1110: }
1114: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1115: {
1116: PetscInt pStart, pEnd, p, localSize = 0;
1120: PetscSectionGetChart(s, &pStart, &pEnd);
1121: for (p = pStart; p < pEnd; ++p) {
1122: PetscInt dof;
1124: PetscSectionGetDof(s, p, &dof);
1125: if (dof > 0) ++localSize;
1126: }
1127: PetscLayoutCreate(comm, layout);
1128: PetscLayoutSetLocalSize(*layout, localSize);
1129: PetscLayoutSetBlockSize(*layout, 1);
1130: PetscLayoutSetUp(*layout);
1131: return(0);
1132: }
1136: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1137: {
1138: PetscInt pStart, pEnd, p, localSize = 0;
1142: PetscSectionGetChart(s, &pStart, &pEnd);
1143: for (p = pStart; p < pEnd; ++p) {
1144: PetscInt dof,cdof;
1146: PetscSectionGetDof(s, p, &dof);
1147: PetscSectionGetConstraintDof(s, p, &cdof);
1148: if (dof-cdof > 0) localSize += dof-cdof;
1149: }
1150: PetscLayoutCreate(comm, layout);
1151: PetscLayoutSetLocalSize(*layout, localSize);
1152: PetscLayoutSetBlockSize(*layout, 1);
1153: PetscLayoutSetUp(*layout);
1154: return(0);
1155: }
1159: /*@
1160: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1162: Not collective
1164: Input Parameters:
1165: + s - the PetscSection
1166: - point - the point
1168: Output Parameter:
1169: . offset - the offset
1171: Level: intermediate
1173: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1174: @*/
1175: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1176: {
1178: #if defined(PETSC_USE_DEBUG)
1179: 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);
1180: #endif
1181: *offset = s->atlasOff[point - s->pStart];
1182: return(0);
1183: }
1187: /*@
1188: PetscSectionSetOffset - Set 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
1195: - offset - the offset
1197: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1199: Level: intermediate
1201: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1202: @*/
1203: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1204: {
1206: 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);
1207: s->atlasOff[point - s->pStart] = offset;
1208: return(0);
1209: }
1213: /*@
1214: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1216: Not collective
1218: Input Parameters:
1219: + s - the PetscSection
1220: . point - the point
1221: - field - the field
1223: Output Parameter:
1224: . offset - the offset
1226: Level: intermediate
1228: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1229: @*/
1230: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1231: {
1235: 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);
1236: PetscSectionGetOffset(s->field[field], point, offset);
1237: return(0);
1238: }
1242: /*@
1243: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1245: Not collective
1247: Input Parameters:
1248: + s - the PetscSection
1249: . point - the point
1250: . field - the field
1251: - offset - the offset
1253: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1255: Level: intermediate
1257: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1258: @*/
1259: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1260: {
1264: 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);
1265: PetscSectionSetOffset(s->field[field], point, offset);
1266: return(0);
1267: }
1271: /* This gives the offset on a point of the field, ignoring constraints */
1272: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1273: {
1274: PetscInt off, foff;
1278: 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);
1279: PetscSectionGetOffset(s, point, &off);
1280: PetscSectionGetOffset(s->field[field], point, &foff);
1281: *offset = foff - off;
1282: return(0);
1283: }
1287: /*@
1288: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1290: Not collective
1292: Input Parameter:
1293: . s - the PetscSection
1295: Output Parameters:
1296: + start - the minimum offset
1297: - end - one more than the maximum offset
1299: Level: intermediate
1301: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1302: @*/
1303: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1304: {
1305: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1309: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1310: PetscSectionGetChart(s, &pStart, &pEnd);
1311: for (p = 0; p < pEnd-pStart; ++p) {
1312: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1314: if (off >= 0) {
1315: os = PetscMin(os, off);
1316: oe = PetscMax(oe, off+dof);
1317: }
1318: }
1319: if (start) *start = os;
1320: if (end) *end = oe;
1321: return(0);
1322: }
1326: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, PetscInt fields[], PetscSection *subs)
1327: {
1328: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1332: if (!numFields) return(0);
1333: PetscSectionGetNumFields(s, &nF);
1334: if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1335: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1336: PetscSectionSetNumFields(*subs, numFields);
1337: for (f = 0; f < numFields; ++f) {
1338: const char *name = NULL;
1339: PetscInt numComp = 0;
1341: PetscSectionGetFieldName(s, fields[f], &name);
1342: PetscSectionSetFieldName(*subs, f, name);
1343: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1344: PetscSectionSetFieldComponents(*subs, f, numComp);
1345: }
1346: PetscSectionGetChart(s, &pStart, &pEnd);
1347: PetscSectionSetChart(*subs, pStart, pEnd);
1348: for (p = pStart; p < pEnd; ++p) {
1349: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1351: for (f = 0; f < numFields; ++f) {
1352: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1353: PetscSectionSetFieldDof(*subs, p, f, fdof);
1354: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1355: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1356: dof += fdof;
1357: cdof += cfdof;
1358: }
1359: PetscSectionSetDof(*subs, p, dof);
1360: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1361: maxCdof = PetscMax(cdof, maxCdof);
1362: }
1363: PetscSectionSetUp(*subs);
1364: if (maxCdof) {
1365: PetscInt *indices;
1367: PetscMalloc1(maxCdof, &indices);
1368: for (p = pStart; p < pEnd; ++p) {
1369: PetscInt cdof;
1371: PetscSectionGetConstraintDof(*subs, p, &cdof);
1372: if (cdof) {
1373: const PetscInt *oldIndices = NULL;
1374: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1376: for (f = 0; f < numFields; ++f) {
1377: PetscInt oldFoff = 0, oldf;
1379: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1380: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1381: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1382: /* This can be sped up if we assume sorted fields */
1383: for (oldf = 0; oldf < fields[f]; ++oldf) {
1384: PetscInt oldfdof = 0;
1385: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1386: oldFoff += oldfdof;
1387: }
1388: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1389: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1390: numConst += cfdof;
1391: fOff += fdof;
1392: }
1393: PetscSectionSetConstraintIndices(*subs, p, indices);
1394: }
1395: }
1396: PetscFree(indices);
1397: }
1398: return(0);
1399: }
1403: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1404: {
1405: const PetscInt *points = NULL, *indices = NULL;
1406: PetscInt numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;
1410: PetscSectionGetNumFields(s, &numFields);
1411: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1412: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1413: for (f = 0; f < numFields; ++f) {
1414: const char *name = NULL;
1415: PetscInt numComp = 0;
1417: PetscSectionGetFieldName(s, f, &name);
1418: PetscSectionSetFieldName(*subs, f, name);
1419: PetscSectionGetFieldComponents(s, f, &numComp);
1420: PetscSectionSetFieldComponents(*subs, f, numComp);
1421: }
1422: /* For right now, we do not try to squeeze the subchart */
1423: if (subpointMap) {
1424: ISGetSize(subpointMap, &numSubpoints);
1425: ISGetIndices(subpointMap, &points);
1426: }
1427: PetscSectionGetChart(s, &pStart, &pEnd);
1428: PetscSectionSetChart(*subs, 0, numSubpoints);
1429: for (p = pStart; p < pEnd; ++p) {
1430: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1432: PetscFindInt(p, numSubpoints, points, &subp);
1433: if (subp < 0) continue;
1434: for (f = 0; f < numFields; ++f) {
1435: PetscSectionGetFieldDof(s, p, f, &fdof);
1436: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1437: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1438: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1439: }
1440: PetscSectionGetDof(s, p, &dof);
1441: PetscSectionSetDof(*subs, subp, dof);
1442: PetscSectionGetConstraintDof(s, p, &cdof);
1443: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1444: }
1445: PetscSectionSetUp(*subs);
1446: /* Change offsets to original offsets */
1447: for (p = pStart; p < pEnd; ++p) {
1448: PetscInt off, foff = 0;
1450: PetscFindInt(p, numSubpoints, points, &subp);
1451: if (subp < 0) continue;
1452: for (f = 0; f < numFields; ++f) {
1453: PetscSectionGetFieldOffset(s, p, f, &foff);
1454: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1455: }
1456: PetscSectionGetOffset(s, p, &off);
1457: PetscSectionSetOffset(*subs, subp, off);
1458: }
1459: /* Copy constraint indices */
1460: for (subp = 0; subp < numSubpoints; ++subp) {
1461: PetscInt cdof;
1463: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1464: if (cdof) {
1465: for (f = 0; f < numFields; ++f) {
1466: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1467: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1468: }
1469: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1470: PetscSectionSetConstraintIndices(*subs, subp, indices);
1471: }
1472: }
1473: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1474: return(0);
1475: }
1479: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1480: {
1481: PetscInt p;
1482: PetscMPIInt rank;
1486: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1487: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
1488: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1489: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1490: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1491: PetscInt b;
1493: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1494: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1495: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1496: }
1497: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1498: } else {
1499: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1500: }
1501: }
1502: PetscViewerFlush(viewer);
1503: return(0);
1504: }
1508: /*@C
1509: PetscSectionView - Views a PetscSection
1511: Collective on PetscSection
1513: Input Parameters:
1514: + s - the PetscSection object to view
1515: - v - the viewer
1517: Level: developer
1519: .seealso PetscSectionCreate(), PetscSectionDestroy()
1520: @*/
1521: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1522: {
1523: PetscBool isascii;
1524: PetscInt f;
1528: if (!viewer) {PetscViewerASCIIGetStdout(PETSC_COMM_SELF, &viewer);}
1530: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1531: if (isascii) {
1532: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1533: if (s->numFields) {
1534: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1535: for (f = 0; f < s->numFields; ++f) {
1536: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
1537: PetscSectionView_ASCII(s->field[f], viewer);
1538: }
1539: } else {
1540: PetscSectionView_ASCII(s, viewer);
1541: }
1542: }
1543: return(0);
1544: }
1548: /*@
1549: PetscSectionReset - Frees all section data.
1551: Not collective
1553: Input Parameters:
1554: . s - the PetscSection
1556: Level: developer
1558: .seealso: PetscSection, PetscSectionCreate()
1559: @*/
1560: PetscErrorCode PetscSectionReset(PetscSection s)
1561: {
1562: PetscInt f;
1566: PetscFree(s->numFieldComponents);
1567: for (f = 0; f < s->numFields; ++f) {
1568: PetscSectionDestroy(&s->field[f]);
1569: PetscFree(s->fieldNames[f]);
1570: }
1571: PetscFree(s->fieldNames);
1572: PetscFree(s->field);
1573: PetscSectionDestroy(&s->bc);
1574: PetscFree(s->bcIndices);
1575: PetscFree2(s->atlasDof, s->atlasOff);
1576: PetscSectionDestroy(&s->clSection);
1577: ISDestroy(&s->clPoints);
1579: s->pStart = -1;
1580: s->pEnd = -1;
1581: s->maxDof = 0;
1582: s->setup = PETSC_FALSE;
1583: s->numFields = 0;
1584: s->clObj = NULL;
1585: return(0);
1586: }
1590: /*@
1591: PetscSectionDestroy - Frees a section object and frees its range if that exists.
1593: Not collective
1595: Input Parameters:
1596: . s - the PetscSection
1598: Level: developer
1600: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1601: recommended they not be used in user codes unless you really gain something in their use.
1603: .seealso: PetscSection, PetscSectionCreate()
1604: @*/
1605: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1606: {
1610: if (!*s) return(0);
1612: if (--((PetscObject)(*s))->refct > 0) {
1613: *s = NULL;
1614: return(0);
1615: }
1616: PetscSectionReset(*s);
1617: PetscHeaderDestroy(s);
1618: return(0);
1619: }
1623: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1624: {
1625: const PetscInt p = point - s->pStart;
1628: *values = &baseArray[s->atlasOff[p]];
1629: return(0);
1630: }
1634: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1635: {
1636: PetscInt *array;
1637: const PetscInt p = point - s->pStart;
1638: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1639: PetscInt cDim = 0;
1643: PetscSectionGetConstraintDof(s, p, &cDim);
1644: array = &baseArray[s->atlasOff[p]];
1645: if (!cDim) {
1646: if (orientation >= 0) {
1647: const PetscInt dim = s->atlasDof[p];
1648: PetscInt i;
1650: if (mode == INSERT_VALUES) {
1651: for (i = 0; i < dim; ++i) array[i] = values[i];
1652: } else {
1653: for (i = 0; i < dim; ++i) array[i] += values[i];
1654: }
1655: } else {
1656: PetscInt offset = 0;
1657: PetscInt j = -1, field, i;
1659: for (field = 0; field < s->numFields; ++field) {
1660: const PetscInt dim = s->field[field]->atlasDof[p];
1662: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1663: offset += dim;
1664: }
1665: }
1666: } else {
1667: if (orientation >= 0) {
1668: const PetscInt dim = s->atlasDof[p];
1669: PetscInt cInd = 0, i;
1670: const PetscInt *cDof;
1672: PetscSectionGetConstraintIndices(s, point, &cDof);
1673: if (mode == INSERT_VALUES) {
1674: for (i = 0; i < dim; ++i) {
1675: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1676: array[i] = values[i];
1677: }
1678: } else {
1679: for (i = 0; i < dim; ++i) {
1680: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1681: array[i] += values[i];
1682: }
1683: }
1684: } else {
1685: const PetscInt *cDof;
1686: PetscInt offset = 0;
1687: PetscInt cOffset = 0;
1688: PetscInt j = 0, field;
1690: PetscSectionGetConstraintIndices(s, point, &cDof);
1691: for (field = 0; field < s->numFields; ++field) {
1692: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
1693: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1694: const PetscInt sDim = dim - tDim;
1695: PetscInt cInd = 0, i,k;
1697: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1698: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1699: array[j] = values[k];
1700: }
1701: offset += dim;
1702: cOffset += dim - tDim;
1703: }
1704: }
1705: }
1706: return(0);
1707: }
1711: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1712: {
1716: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1717: return(0);
1718: }
1722: /*@C
1723: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
1725: Input Parameters:
1726: + s - The PetscSection
1727: - point - The point
1729: Output Parameter:
1730: . indices - The constrained dofs
1732: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
1734: Level: advanced
1736: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1737: @*/
1738: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1739: {
1743: if (s->bc) {
1744: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1745: } else *indices = NULL;
1746: return(0);
1747: }
1751: /*@C
1752: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
1754: Input Parameters:
1755: + s - The PetscSection
1756: . point - The point
1757: - indices - The constrained dofs
1759: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
1761: Level: advanced
1763: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1764: @*/
1765: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1766: {
1770: if (s->bc) {
1771: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1772: }
1773: return(0);
1774: }
1778: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1779: {
1783: 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);
1784: PetscSectionGetConstraintIndices(s->field[field], point, indices);
1785: return(0);
1786: }
1790: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
1791: {
1795: 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);
1796: PetscSectionSetConstraintIndices(s->field[field], point, indices);
1797: return(0);
1798: }
1802: /*@
1803: PetscSectionPermute - Reorder the section according to the input point permutation
1805: Collective on PetscSection
1807: Input Parameter:
1808: + section - The PetscSection object
1809: - perm - The point permutation, old point p becomes new point perm[p]
1811: Output Parameter:
1812: . sectionNew - The permuted PetscSection
1814: Level: intermediate
1816: .keywords: mesh
1817: .seealso: MatPermute()
1818: @*/
1819: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1820: {
1821: PetscSection s = section, sNew;
1822: const PetscInt *perm;
1823: PetscInt numFields, f, numPoints, pStart, pEnd, p;
1824: PetscErrorCode ierr;
1830: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1831: PetscSectionGetNumFields(s, &numFields);
1832: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1833: for (f = 0; f < numFields; ++f) {
1834: const char *name;
1835: PetscInt numComp;
1837: PetscSectionGetFieldName(s, f, &name);
1838: PetscSectionSetFieldName(sNew, f, name);
1839: PetscSectionGetFieldComponents(s, f, &numComp);
1840: PetscSectionSetFieldComponents(sNew, f, numComp);
1841: }
1842: ISGetLocalSize(permutation, &numPoints);
1843: ISGetIndices(permutation, &perm);
1844: PetscSectionGetChart(s, &pStart, &pEnd);
1845: PetscSectionSetChart(sNew, pStart, pEnd);
1846: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1847: for (p = pStart; p < pEnd; ++p) {
1848: PetscInt dof, cdof;
1850: PetscSectionGetDof(s, p, &dof);
1851: PetscSectionSetDof(sNew, perm[p], dof);
1852: PetscSectionGetConstraintDof(s, p, &cdof);
1853: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1854: for (f = 0; f < numFields; ++f) {
1855: PetscSectionGetFieldDof(s, p, f, &dof);
1856: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1857: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1858: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1859: }
1860: }
1861: PetscSectionSetUp(sNew);
1862: for (p = pStart; p < pEnd; ++p) {
1863: const PetscInt *cind;
1864: PetscInt cdof;
1866: PetscSectionGetConstraintDof(s, p, &cdof);
1867: if (cdof) {
1868: PetscSectionGetConstraintIndices(s, p, &cind);
1869: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1870: }
1871: for (f = 0; f < numFields; ++f) {
1872: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1873: if (cdof) {
1874: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1875: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1876: }
1877: }
1878: }
1879: ISRestoreIndices(permutation, &perm);
1880: *sectionNew = sNew;
1881: return(0);
1882: }
1884: /*
1885: I need extract and merge routines for section based on fields. This should be trivial except for updating the
1886: constraint indices.
1888: Then I need a new interface for DMCreateDecomposition that takes groups of fields and returns a real DMPlex
1889: that shares the mesh parts, and has the extracted section
1890: */
1894: PetscErrorCode PetscSFConvertPartition(PetscSF sfPart, PetscSection partSection, IS partition, ISLocalToGlobalMapping *renumbering, PetscSF *sf)
1895: {
1896: MPI_Comm comm;
1897: PetscSF sfPoints;
1898: PetscInt *partSizes,*partOffsets,p,i,numParts,numMyPoints,numPoints,count;
1899: const PetscInt *partArray;
1900: PetscSFNode *sendPoints;
1901: PetscMPIInt rank;
1905: PetscObjectGetComm((PetscObject)sfPart,&comm);
1906: MPI_Comm_rank(comm, &rank);
1908: /* Get the number of parts and sizes that I have to distribute */
1909: PetscSFGetGraph(sfPart,NULL,&numParts,NULL,NULL);
1910: PetscMalloc2(numParts,&partSizes,numParts,&partOffsets);
1911: for (p=0,numPoints=0; p<numParts; p++) {
1912: PetscSectionGetDof(partSection, p, &partSizes[p]);
1913: numPoints += partSizes[p];
1914: }
1915: numMyPoints = 0;
1916: PetscSFFetchAndOpBegin(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1917: PetscSFFetchAndOpEnd(sfPart,MPIU_INT,&numMyPoints,partSizes,partOffsets,MPIU_SUM);
1918: /* I will receive numMyPoints. I will send a total of numPoints, to be placed on remote procs at partOffsets */
1920: /* Create SF mapping locations (addressed through partition, as indexed leaves) to new owners (roots) */
1921: PetscMalloc1(numPoints,&sendPoints);
1922: for (p=0,count=0; p<numParts; p++) {
1923: for (i=0; i<partSizes[p]; i++) {
1924: sendPoints[count].rank = p;
1925: sendPoints[count].index = partOffsets[p]+i;
1926: count++;
1927: }
1928: }
1929: if (count != numPoints) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count %D should equal numPoints=%D",count,numPoints);
1930: PetscFree2(partSizes,partOffsets);
1931: ISGetIndices(partition,&partArray);
1932: PetscSFCreate(comm,&sfPoints);
1933: PetscSFSetFromOptions(sfPoints);
1934: PetscSFSetGraph(sfPoints,numMyPoints,numPoints,partArray,PETSC_USE_POINTER,sendPoints,PETSC_OWN_POINTER);
1936: /* Invert SF so that the new owners are leaves and the locations indexed through partition are the roots */
1937: PetscSFCreateInverseSF(sfPoints,sf);
1938: PetscSFDestroy(&sfPoints);
1939: ISRestoreIndices(partition,&partArray);
1941: /* Create the new local-to-global mapping */
1942: ISLocalToGlobalMappingCreateSF(*sf,0,renumbering);
1943: return(0);
1944: }
1948: /*@C
1949: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
1951: Collective
1953: Input Parameters:
1954: + sf - The SF
1955: - rootSection - Section defined on root space
1957: Output Parameters:
1958: + remoteOffsets - root offsets in leaf storage, or NULL
1959: - leafSection - Section defined on the leaf space
1961: Level: intermediate
1963: .seealso: PetscSFCreate()
1964: @*/
1965: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1966: {
1967: PetscSF embedSF;
1968: const PetscInt *ilocal, *indices;
1969: IS selected;
1970: PetscInt numFields, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;
1974: PetscSectionGetNumFields(rootSection, &numFields);
1975: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
1976: for (f = 0; f < numFields; ++f) {
1977: PetscInt numComp = 0;
1978: PetscSectionGetFieldComponents(rootSection, f, &numComp);
1979: PetscSectionSetFieldComponents(leafSection, f, numComp);
1980: }
1981: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1982: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1983: ISGetIndices(selected, &indices);
1984: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1985: ISRestoreIndices(selected, &indices);
1986: ISDestroy(&selected);
1987: PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
1988: if (nleaves && ilocal) {
1989: for (i = 0; i < nleaves; ++i) {
1990: lpStart = PetscMin(lpStart, ilocal[i]);
1991: lpEnd = PetscMax(lpEnd, ilocal[i]);
1992: }
1993: ++lpEnd;
1994: } else {
1995: lpStart = 0;
1996: lpEnd = nleaves;
1997: }
1998: PetscSectionSetChart(leafSection, lpStart, lpEnd);
1999: /* Could fuse these at the cost of a copy and extra allocation */
2000: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2001: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2002: for (f = 0; f < numFields; ++f) {
2003: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2004: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2005: }
2006: if (remoteOffsets) {
2007: PetscMalloc1((lpEnd - lpStart), remoteOffsets);
2008: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2009: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2010: }
2011: PetscSFDestroy(&embedSF);
2012: PetscSectionSetUp(leafSection);
2013: return(0);
2014: }
2018: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2019: {
2020: PetscSF embedSF;
2021: const PetscInt *indices;
2022: IS selected;
2023: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2024: PetscErrorCode ierr;
2027: *remoteOffsets = NULL;
2028: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2029: if (numRoots < 0) return(0);
2030: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2031: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2032: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2033: ISGetIndices(selected, &indices);
2034: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2035: ISRestoreIndices(selected, &indices);
2036: ISDestroy(&selected);
2037: PetscCalloc1((lpEnd - lpStart), remoteOffsets);
2038: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2039: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2040: PetscSFDestroy(&embedSF);
2041: return(0);
2042: }
2046: /*@C
2047: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2049: Input Parameters:
2050: + sf - The SF
2051: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section), or NULL
2052: - remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2054: Output Parameters:
2055: + leafSection - Data layout of local points for incoming data (this is the distributed section)
2056: - sectionSF - The new SF
2058: Note: Either rootSection or remoteOffsets can be specified
2060: Level: intermediate
2062: .seealso: PetscSFCreate()
2063: @*/
2064: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2065: {
2066: MPI_Comm comm;
2067: const PetscInt *localPoints;
2068: const PetscSFNode *remotePoints;
2069: PetscInt lpStart, lpEnd;
2070: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2071: PetscInt *localIndices;
2072: PetscSFNode *remoteIndices;
2073: PetscInt i, ind;
2074: PetscErrorCode ierr;
2082: PetscObjectGetComm((PetscObject)sf,&comm);
2083: PetscSFCreate(comm, sectionSF);
2084: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2085: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2086: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2087: if (numRoots < 0) return(0);
2088: for (i = 0; i < numPoints; ++i) {
2089: PetscInt localPoint = localPoints ? localPoints[i] : i;
2090: PetscInt dof;
2092: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2093: PetscSectionGetDof(leafSection, localPoint, &dof);
2094: numIndices += dof;
2095: }
2096: }
2097: PetscMalloc1(numIndices, &localIndices);
2098: PetscMalloc1(numIndices, &remoteIndices);
2099: /* Create new index graph */
2100: for (i = 0, ind = 0; i < numPoints; ++i) {
2101: PetscInt localPoint = localPoints ? localPoints[i] : i;
2102: PetscInt rank = remotePoints[i].rank;
2104: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2105: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2106: PetscInt loff, dof, d;
2108: PetscSectionGetOffset(leafSection, localPoint, &loff);
2109: PetscSectionGetDof(leafSection, localPoint, &dof);
2110: for (d = 0; d < dof; ++d, ++ind) {
2111: localIndices[ind] = loff+d;
2112: remoteIndices[ind].rank = rank;
2113: remoteIndices[ind].index = remoteOffset+d;
2114: }
2115: }
2116: }
2117: PetscFree(remoteOffsets);
2118: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2119: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2120: return(0);
2121: }
2125: /*@
2126: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2128: Input Parameters:
2129: + section - The PetscSection
2130: . obj - A PetscObject which serves as the key for this index
2131: . clSection - Section giving the size of the closure of each point
2132: - clPoints - IS giving the points in each closure
2134: Note: We compress out closure points with no dofs in this section
2136: Level: intermediate
2138: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2139: @*/
2140: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2141: {
2145: section->clObj = obj;
2146: PetscSectionDestroy(§ion->clSection);
2147: ISDestroy(§ion->clPoints);
2148: section->clSection = clSection;
2149: section->clPoints = clPoints;
2150: return(0);
2151: }
2155: /*@
2156: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2158: Input Parameters:
2159: + section - The PetscSection
2160: - obj - A PetscObject which serves as the key for this index
2162: Output Parameters:
2163: + clSection - Section giving the size of the closure of each point
2164: - clPoints - IS giving the points in each closure
2166: Note: We compress out closure points with no dofs in this section
2168: Level: intermediate
2170: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2171: @*/
2172: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2173: {
2175: if (section->clObj == obj) {
2176: if (clSection) *clSection = section->clSection;
2177: if (clPoints) *clPoints = section->clPoints;
2178: } else {
2179: if (clSection) *clSection = NULL;
2180: if (clPoints) *clPoints = NULL;
2181: }
2182: return(0);
2183: }
2187: /*@
2188: PetscSectionGetField - Get the subsection associated with a single field
2190: Input Parameters:
2191: + s - The PetscSection
2192: - field - The field number
2194: Output Parameter:
2195: . subs - The subsection for the given field
2197: Level: intermediate
2199: .seealso: PetscSectionSetNumFields()
2200: @*/
2201: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2202: {
2208: 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);
2209: PetscObjectReference((PetscObject) s->field[field]);
2210: *subs = s->field[field];
2211: return(0);
2212: }