Actual source code: section.c
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsf.h>
8: PetscClassId PETSC_SECTION_CLASSID;
10: /*@
11: PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
13: Collective
15: Input Parameters:
16: + comm - the MPI communicator
17: - s - pointer to the section
19: Level: beginner
21: Notes:
22: Typical calling sequence
23: $ PetscSectionCreate(MPI_Comm,PetscSection *);
24: $ PetscSectionSetNumFields(PetscSection, numFields);
25: $ PetscSectionSetChart(PetscSection,low,high);
26: $ PetscSectionSetDof(PetscSection,point,numdof);
27: $ PetscSectionSetUp(PetscSection);
28: $ PetscSectionGetOffset(PetscSection,point,PetscInt *);
29: $ PetscSectionDestroy(PetscSection);
31: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementations; it is
32: recommended they not be used in user codes unless you really gain something in their use.
34: .seealso: PetscSection, PetscSectionDestroy()
35: @*/
36: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
37: {
39: ISInitializePackage();
41: PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
43: (*s)->pStart = -1;
44: (*s)->pEnd = -1;
45: (*s)->perm = NULL;
46: (*s)->pointMajor = PETSC_TRUE;
47: (*s)->includesConstraints = PETSC_TRUE;
48: (*s)->atlasDof = NULL;
49: (*s)->atlasOff = NULL;
50: (*s)->bc = NULL;
51: (*s)->bcIndices = NULL;
52: (*s)->setup = PETSC_FALSE;
53: (*s)->numFields = 0;
54: (*s)->fieldNames = NULL;
55: (*s)->field = NULL;
56: (*s)->useFieldOff = PETSC_FALSE;
57: (*s)->compNames = NULL;
58: (*s)->clObj = NULL;
59: (*s)->clHash = NULL;
60: (*s)->clSection = NULL;
61: (*s)->clPoints = NULL;
62: PetscSectionInvalidateMaxDof_Internal(*s);
63: return 0;
64: }
66: /*@
67: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
69: Collective
71: Input Parameter:
72: . section - the PetscSection
74: Output Parameter:
75: . newSection - the copy
77: Level: intermediate
79: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
80: @*/
81: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
82: {
83: PetscSectionSym sym;
84: IS perm;
85: PetscInt numFields, f, c, pStart, pEnd, p;
89: PetscSectionReset(newSection);
90: PetscSectionGetNumFields(section, &numFields);
91: if (numFields) PetscSectionSetNumFields(newSection, numFields);
92: for (f = 0; f < numFields; ++f) {
93: const char *fieldName = NULL, *compName = NULL;
94: PetscInt numComp = 0;
96: PetscSectionGetFieldName(section, f, &fieldName);
97: PetscSectionSetFieldName(newSection, f, fieldName);
98: PetscSectionGetFieldComponents(section, f, &numComp);
99: PetscSectionSetFieldComponents(newSection, f, numComp);
100: for (c = 0; c < numComp; ++c) {
101: PetscSectionGetComponentName(section, f, c, &compName);
102: PetscSectionSetComponentName(newSection, f, c, compName);
103: }
104: PetscSectionGetFieldSym(section, f, &sym);
105: PetscSectionSetFieldSym(newSection, f, sym);
106: }
107: PetscSectionGetChart(section, &pStart, &pEnd);
108: PetscSectionSetChart(newSection, pStart, pEnd);
109: PetscSectionGetPermutation(section, &perm);
110: PetscSectionSetPermutation(newSection, perm);
111: PetscSectionGetSym(section, &sym);
112: PetscSectionSetSym(newSection, sym);
113: for (p = pStart; p < pEnd; ++p) {
114: PetscInt dof, cdof, fcdof = 0;
116: PetscSectionGetDof(section, p, &dof);
117: PetscSectionSetDof(newSection, p, dof);
118: PetscSectionGetConstraintDof(section, p, &cdof);
119: if (cdof) PetscSectionSetConstraintDof(newSection, p, cdof);
120: for (f = 0; f < numFields; ++f) {
121: PetscSectionGetFieldDof(section, p, f, &dof);
122: PetscSectionSetFieldDof(newSection, p, f, dof);
123: if (cdof) {
124: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
125: if (fcdof) PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);
126: }
127: }
128: }
129: PetscSectionSetUp(newSection);
130: for (p = pStart; p < pEnd; ++p) {
131: PetscInt off, cdof, fcdof = 0;
132: const PetscInt *cInd;
134: /* Must set offsets in case they do not agree with the prefix sums */
135: PetscSectionGetOffset(section, p, &off);
136: PetscSectionSetOffset(newSection, p, off);
137: PetscSectionGetConstraintDof(section, p, &cdof);
138: if (cdof) {
139: PetscSectionGetConstraintIndices(section, p, &cInd);
140: PetscSectionSetConstraintIndices(newSection, p, cInd);
141: for (f = 0; f < numFields; ++f) {
142: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
143: if (fcdof) {
144: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
145: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
146: }
147: }
148: }
149: }
150: return 0;
151: }
153: /*@
154: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
156: Collective
158: Input Parameter:
159: . section - the PetscSection
161: Output Parameter:
162: . newSection - the copy
164: Level: beginner
166: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
167: @*/
168: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
169: {
172: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
173: PetscSectionCopy(section, *newSection);
174: return 0;
175: }
177: /*@
178: PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database
180: Collective
182: Input Parameter:
183: . section - the PetscSection
185: Options Database:
186: . -petscsection_point_major - PETSC_TRUE for point-major order
188: Level: intermediate
190: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
191: @*/
192: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
193: {
197: PetscObjectOptionsBegin((PetscObject) s);
198: PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
199: /* process any options handlers added with PetscObjectAddOptionsHandler() */
200: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
201: PetscOptionsEnd();
202: PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
203: return 0;
204: }
206: /*@
207: PetscSectionCompare - Compares two sections
209: Collective
211: Input Parameters:
212: + s1 - the first PetscSection
213: - s2 - the second PetscSection
215: Output Parameter:
216: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
218: Level: intermediate
220: Notes:
221: Field names are disregarded.
223: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
224: @*/
225: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
226: {
227: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
228: const PetscInt *idx1, *idx2;
229: IS perm1, perm2;
230: PetscBool flg;
231: PetscMPIInt mflg;
236: flg = PETSC_FALSE;
238: MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
239: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
240: *congruent = PETSC_FALSE;
241: return 0;
242: }
244: PetscSectionGetChart(s1, &pStart, &pEnd);
245: PetscSectionGetChart(s2, &n1, &n2);
246: if (pStart != n1 || pEnd != n2) goto not_congruent;
248: PetscSectionGetPermutation(s1, &perm1);
249: PetscSectionGetPermutation(s2, &perm2);
250: if (perm1 && perm2) {
251: ISEqual(perm1, perm2, congruent);
252: if (!(*congruent)) goto not_congruent;
253: } else if (perm1 != perm2) goto not_congruent;
255: for (p = pStart; p < pEnd; ++p) {
256: PetscSectionGetOffset(s1, p, &n1);
257: PetscSectionGetOffset(s2, p, &n2);
258: if (n1 != n2) goto not_congruent;
260: PetscSectionGetDof(s1, p, &n1);
261: PetscSectionGetDof(s2, p, &n2);
262: if (n1 != n2) goto not_congruent;
264: PetscSectionGetConstraintDof(s1, p, &ncdof);
265: PetscSectionGetConstraintDof(s2, p, &n2);
266: if (ncdof != n2) goto not_congruent;
268: PetscSectionGetConstraintIndices(s1, p, &idx1);
269: PetscSectionGetConstraintIndices(s2, p, &idx2);
270: PetscArraycmp(idx1, idx2, ncdof, congruent);
271: if (!(*congruent)) goto not_congruent;
272: }
274: PetscSectionGetNumFields(s1, &nfields);
275: PetscSectionGetNumFields(s2, &n2);
276: if (nfields != n2) goto not_congruent;
278: for (f = 0; f < nfields; ++f) {
279: PetscSectionGetFieldComponents(s1, f, &n1);
280: PetscSectionGetFieldComponents(s2, f, &n2);
281: if (n1 != n2) goto not_congruent;
283: for (p = pStart; p < pEnd; ++p) {
284: PetscSectionGetFieldOffset(s1, p, f, &n1);
285: PetscSectionGetFieldOffset(s2, p, f, &n2);
286: if (n1 != n2) goto not_congruent;
288: PetscSectionGetFieldDof(s1, p, f, &n1);
289: PetscSectionGetFieldDof(s2, p, f, &n2);
290: if (n1 != n2) goto not_congruent;
292: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
293: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
294: if (nfcdof != n2) goto not_congruent;
296: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
297: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
298: PetscArraycmp(idx1, idx2, nfcdof, congruent);
299: if (!(*congruent)) goto not_congruent;
300: }
301: }
303: flg = PETSC_TRUE;
304: not_congruent:
305: MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
306: return 0;
307: }
309: /*@
310: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
312: Not Collective
314: Input Parameter:
315: . s - the PetscSection
317: Output Parameter:
318: . numFields - the number of fields defined, or 0 if none were defined
320: Level: intermediate
322: .seealso: PetscSectionSetNumFields()
323: @*/
324: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
325: {
328: *numFields = s->numFields;
329: return 0;
330: }
332: /*@
333: PetscSectionSetNumFields - Sets the number of fields.
335: Not Collective
337: Input Parameters:
338: + s - the PetscSection
339: - numFields - the number of fields
341: Level: intermediate
343: .seealso: PetscSectionGetNumFields()
344: @*/
345: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
346: {
347: PetscInt f;
351: PetscSectionReset(s);
353: s->numFields = numFields;
354: PetscMalloc1(s->numFields, &s->numFieldComponents);
355: PetscMalloc1(s->numFields, &s->fieldNames);
356: PetscMalloc1(s->numFields, &s->compNames);
357: PetscMalloc1(s->numFields, &s->field);
358: for (f = 0; f < s->numFields; ++f) {
359: char name[64];
361: s->numFieldComponents[f] = 1;
363: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
364: PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f);
365: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
366: PetscSNPrintf(name, 64, "Component_0");
367: PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
368: PetscStrallocpy(name, (char **) &s->compNames[f][0]);
369: }
370: return 0;
371: }
373: /*@C
374: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
376: Not Collective
378: Input Parameters:
379: + s - the PetscSection
380: - field - the field number
382: Output Parameter:
383: . fieldName - the field name
385: Level: intermediate
387: .seealso: PetscSectionSetFieldName()
388: @*/
389: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
390: {
393: PetscSectionCheckValidField(field,s->numFields);
394: *fieldName = s->fieldNames[field];
395: return 0;
396: }
398: /*@C
399: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
401: Not Collective
403: Input Parameters:
404: + s - the PetscSection
405: . field - the field number
406: - fieldName - the field name
408: Level: intermediate
410: .seealso: PetscSectionGetFieldName()
411: @*/
412: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
413: {
416: PetscSectionCheckValidField(field,s->numFields);
417: PetscFree(s->fieldNames[field]);
418: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
419: return 0;
420: }
422: /*@C
423: PetscSectionGetComponentName - Gets the name of a field component in the PetscSection
425: Not Collective
427: Input Parameters:
428: + s - the PetscSection
429: . field - the field number
430: . comp - the component number
431: - compName - the component name
433: Level: intermediate
435: .seealso: PetscSectionSetComponentName()
436: @*/
437: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
438: {
441: PetscSectionCheckValidField(field,s->numFields);
442: PetscSectionCheckValidFieldComponent(comp,s->numFieldComponents[field]);
443: *compName = s->compNames[field][comp];
444: return 0;
445: }
447: /*@C
448: PetscSectionSetComponentName - Sets the name of a field component in the PetscSection
450: Not Collective
452: Input Parameters:
453: + s - the PetscSection
454: . field - the field number
455: . comp - the component number
456: - compName - the component name
458: Level: intermediate
460: .seealso: PetscSectionGetComponentName()
461: @*/
462: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
463: {
466: PetscSectionCheckValidField(field,s->numFields);
467: PetscSectionCheckValidFieldComponent(comp,s->numFieldComponents[field]);
468: PetscFree(s->compNames[field][comp]);
469: PetscStrallocpy(compName, (char**) &s->compNames[field][comp]);
470: return 0;
471: }
473: /*@
474: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
476: Not Collective
478: Input Parameters:
479: + s - the PetscSection
480: - field - the field number
482: Output Parameter:
483: . numComp - the number of field components
485: Level: intermediate
487: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
488: @*/
489: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
490: {
493: PetscSectionCheckValidField(field,s->numFields);
494: *numComp = s->numFieldComponents[field];
495: return 0;
496: }
498: /*@
499: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
501: Not Collective
503: Input Parameters:
504: + s - the PetscSection
505: . field - the field number
506: - numComp - the number of field components
508: Level: intermediate
510: .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
511: @*/
512: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
513: {
514: PetscInt c;
517: PetscSectionCheckValidField(field,s->numFields);
518: if (s->compNames) {
519: for (c = 0; c < s->numFieldComponents[field]; ++c) {
520: PetscFree(s->compNames[field][c]);
521: }
522: PetscFree(s->compNames[field]);
523: }
525: s->numFieldComponents[field] = numComp;
526: if (numComp) {
527: PetscMalloc1(numComp, (char ***) &s->compNames[field]);
528: for (c = 0; c < numComp; ++c) {
529: char name[64];
531: PetscSNPrintf(name, 64, "%" PetscInt_FMT, c);
532: PetscStrallocpy(name, (char **) &s->compNames[field][c]);
533: }
534: }
535: return 0;
536: }
538: /*@
539: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
541: Not Collective
543: Input Parameter:
544: . s - the PetscSection
546: Output Parameters:
547: + pStart - the first point
548: - pEnd - one past the last point
550: Level: intermediate
552: .seealso: PetscSectionSetChart(), PetscSectionCreate()
553: @*/
554: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
555: {
557: if (pStart) *pStart = s->pStart;
558: if (pEnd) *pEnd = s->pEnd;
559: return 0;
560: }
562: /*@
563: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
565: Not Collective
567: Input Parameters:
568: + s - the PetscSection
569: . pStart - the first point
570: - pEnd - one past the last point
572: Level: intermediate
574: .seealso: PetscSectionGetChart(), PetscSectionCreate()
575: @*/
576: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
577: {
578: PetscInt f;
581: if (pStart == s->pStart && pEnd == s->pEnd) return 0;
582: /* Cannot Reset() because it destroys field information */
583: s->setup = PETSC_FALSE;
584: PetscSectionDestroy(&s->bc);
585: PetscFree(s->bcIndices);
586: PetscFree2(s->atlasDof, s->atlasOff);
588: s->pStart = pStart;
589: s->pEnd = pEnd;
590: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
591: PetscArrayzero(s->atlasDof, pEnd - pStart);
592: for (f = 0; f < s->numFields; ++f) {
593: PetscSectionSetChart(s->field[f], pStart, pEnd);
594: }
595: return 0;
596: }
598: /*@
599: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
601: Not Collective
603: Input Parameter:
604: . s - the PetscSection
606: Output Parameters:
607: . perm - The permutation as an IS
609: Level: intermediate
611: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
612: @*/
613: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
614: {
617: return 0;
618: }
620: /*@
621: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
623: Not Collective
625: Input Parameters:
626: + s - the PetscSection
627: - perm - the permutation of points
629: Level: intermediate
631: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
632: @*/
633: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
634: {
638: if (s->perm != perm) {
639: ISDestroy(&s->perm);
640: if (perm) {
641: s->perm = perm;
642: PetscObjectReference((PetscObject) s->perm);
643: }
644: }
645: return 0;
646: }
648: /*@
649: PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major
651: Not Collective
653: Input Parameter:
654: . s - the PetscSection
656: Output Parameter:
657: . pm - the flag for point major ordering
659: Level: intermediate
661: .seealso: PetscSectionSetPointMajor()
662: @*/
663: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
664: {
667: *pm = s->pointMajor;
668: return 0;
669: }
671: /*@
672: PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major
674: Not Collective
676: Input Parameters:
677: + s - the PetscSection
678: - pm - the flag for point major ordering
680: Level: intermediate
682: .seealso: PetscSectionGetPointMajor()
683: @*/
684: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
685: {
688: s->pointMajor = pm;
689: return 0;
690: }
692: /*@
693: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets
695: Not Collective
697: Input Parameter:
698: . s - the PetscSection
700: Output Parameter:
701: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
703: Level: intermediate
705: .seealso: PetscSectionSetIncludesConstraints()
706: @*/
707: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
708: {
711: *includesConstraints = s->includesConstraints;
712: return 0;
713: }
715: /*@
716: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
718: Not Collective
720: Input Parameters:
721: + s - the PetscSection
722: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
724: Level: intermediate
726: .seealso: PetscSectionGetIncludesConstraints()
727: @*/
728: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
729: {
732: s->includesConstraints = includesConstraints;
733: return 0;
734: }
736: /*@
737: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
739: Not Collective
741: Input Parameters:
742: + s - the PetscSection
743: - point - the point
745: Output Parameter:
746: . numDof - the number of dof
748: Level: intermediate
750: .seealso: PetscSectionSetDof(), PetscSectionCreate()
751: @*/
752: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
753: {
757: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
758: *numDof = s->atlasDof[point - s->pStart];
759: return 0;
760: }
762: /*@
763: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
765: Not Collective
767: Input Parameters:
768: + s - the PetscSection
769: . point - the point
770: - numDof - the number of dof
772: Level: intermediate
774: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
775: @*/
776: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
777: {
779: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
780: s->atlasDof[point - s->pStart] = numDof;
781: PetscSectionInvalidateMaxDof_Internal(s);
782: return 0;
783: }
785: /*@
786: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
788: Not Collective
790: Input Parameters:
791: + s - the PetscSection
792: . point - the point
793: - numDof - the number of additional dof
795: Level: intermediate
797: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
798: @*/
799: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
800: {
803: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
804: s->atlasDof[point - s->pStart] += numDof;
805: PetscSectionInvalidateMaxDof_Internal(s);
806: return 0;
807: }
809: /*@
810: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
812: Not Collective
814: Input Parameters:
815: + s - the PetscSection
816: . point - the point
817: - field - the field
819: Output Parameter:
820: . numDof - the number of dof
822: Level: intermediate
824: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
825: @*/
826: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
827: {
830: PetscSectionCheckValidField(field,s->numFields);
831: PetscSectionGetDof(s->field[field], point, numDof);
832: return 0;
833: }
835: /*@
836: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
838: Not Collective
840: Input Parameters:
841: + s - the PetscSection
842: . point - the point
843: . field - the field
844: - numDof - the number of dof
846: Level: intermediate
848: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
849: @*/
850: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
851: {
853: PetscSectionCheckValidField(field,s->numFields);
854: PetscSectionSetDof(s->field[field], point, numDof);
855: return 0;
856: }
858: /*@
859: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
861: Not Collective
863: Input Parameters:
864: + s - the PetscSection
865: . point - the point
866: . field - the field
867: - numDof - the number of dof
869: Level: intermediate
871: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
872: @*/
873: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
874: {
876: PetscSectionCheckValidField(field,s->numFields);
877: PetscSectionAddDof(s->field[field], point, numDof);
878: return 0;
879: }
881: /*@
882: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
884: Not Collective
886: Input Parameters:
887: + s - the PetscSection
888: - point - the point
890: Output Parameter:
891: . numDof - the number of dof which are fixed by constraints
893: Level: intermediate
895: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
896: @*/
897: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
898: {
901: if (s->bc) {
902: PetscSectionGetDof(s->bc, point, numDof);
903: } else *numDof = 0;
904: return 0;
905: }
907: /*@
908: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
910: Not Collective
912: Input Parameters:
913: + s - the PetscSection
914: . point - the point
915: - numDof - the number of dof which are fixed by constraints
917: Level: intermediate
919: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
920: @*/
921: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
922: {
924: if (numDof) {
925: PetscSectionCheckConstraints_Private(s);
926: PetscSectionSetDof(s->bc, point, numDof);
927: }
928: return 0;
929: }
931: /*@
932: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
934: Not Collective
936: Input Parameters:
937: + s - the PetscSection
938: . point - the point
939: - numDof - the number of additional dof which are fixed by constraints
941: Level: intermediate
943: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
944: @*/
945: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
946: {
948: if (numDof) {
949: PetscSectionCheckConstraints_Private(s);
950: PetscSectionAddDof(s->bc, point, numDof);
951: }
952: return 0;
953: }
955: /*@
956: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
958: Not Collective
960: Input Parameters:
961: + s - the PetscSection
962: . point - the point
963: - field - the field
965: Output Parameter:
966: . numDof - the number of dof which are fixed by constraints
968: Level: intermediate
970: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
971: @*/
972: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
973: {
976: PetscSectionCheckValidField(field,s->numFields);
977: PetscSectionGetConstraintDof(s->field[field], point, numDof);
978: return 0;
979: }
981: /*@
982: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
984: Not Collective
986: Input Parameters:
987: + s - the PetscSection
988: . point - the point
989: . field - the field
990: - numDof - the number of dof which are fixed by constraints
992: Level: intermediate
994: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
995: @*/
996: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
997: {
999: PetscSectionCheckValidField(field,s->numFields);
1000: PetscSectionSetConstraintDof(s->field[field], point, numDof);
1001: return 0;
1002: }
1004: /*@
1005: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1007: Not Collective
1009: Input Parameters:
1010: + s - the PetscSection
1011: . point - the point
1012: . field - the field
1013: - numDof - the number of additional dof which are fixed by constraints
1015: Level: intermediate
1017: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1018: @*/
1019: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1020: {
1022: PetscSectionCheckValidField(field,s->numFields);
1023: PetscSectionAddConstraintDof(s->field[field], point, numDof);
1024: return 0;
1025: }
1027: /*@
1028: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1030: Not Collective
1032: Input Parameter:
1033: . s - the PetscSection
1035: Level: advanced
1037: .seealso: PetscSectionSetUp(), PetscSectionCreate()
1038: @*/
1039: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1040: {
1042: if (s->bc) {
1043: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
1045: PetscSectionSetUp(s->bc);
1046: PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1047: }
1048: return 0;
1049: }
1051: /*@
1052: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
1054: Not Collective
1056: Input Parameter:
1057: . s - the PetscSection
1059: Level: intermediate
1061: .seealso: PetscSectionCreate()
1062: @*/
1063: PetscErrorCode PetscSectionSetUp(PetscSection s)
1064: {
1065: const PetscInt *pind = NULL;
1066: PetscInt offset = 0, foff, p, f;
1069: if (s->setup) return 0;
1070: s->setup = PETSC_TRUE;
1071: /* Set offsets and field offsets for all points */
1072: /* Assume that all fields have the same chart */
1074: if (s->perm) ISGetIndices(s->perm, &pind);
1075: if (s->pointMajor) {
1076: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1077: const PetscInt q = pind ? pind[p] : p;
1079: /* Set point offset */
1080: s->atlasOff[q] = offset;
1081: offset += s->atlasDof[q];
1082: /* Set field offset */
1083: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1084: PetscSection sf = s->field[f];
1086: sf->atlasOff[q] = foff;
1087: foff += sf->atlasDof[q];
1088: }
1089: }
1090: } else {
1091: /* Set field offsets for all points */
1092: for (f = 0; f < s->numFields; ++f) {
1093: PetscSection sf = s->field[f];
1095: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1096: const PetscInt q = pind ? pind[p] : p;
1098: sf->atlasOff[q] = offset;
1099: offset += sf->atlasDof[q];
1100: }
1101: }
1102: /* Disable point offsets since these are unused */
1103: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1104: s->atlasOff[p] = -1;
1105: }
1106: }
1107: if (s->perm) ISRestoreIndices(s->perm, &pind);
1108: /* Setup BC sections */
1109: PetscSectionSetUpBC(s);
1110: for (f = 0; f < s->numFields; ++f) PetscSectionSetUpBC(s->field[f]);
1111: return 0;
1112: }
1114: /*@
1115: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1117: Not Collective
1119: Input Parameters:
1120: . s - the PetscSection
1122: Output Parameter:
1123: . maxDof - the maximum dof
1125: Level: intermediate
1127: Notes:
1128: The returned number is up-to-date without need for PetscSectionSetUp().
1130: Developer Notes:
1131: The returned number is calculated lazily and stashed.
1132: A call to PetscSectionInvalidateMaxDof_Internal() invalidates the stashed value.
1133: PetscSectionInvalidateMaxDof_Internal() is called in PetscSectionSetDof(), PetscSectionAddDof() and PetscSectionReset().
1134: It should also be called every time atlasDof is modified directly.
1136: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionAddDof(), PetscSectionCreate()
1137: @*/
1138: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1139: {
1140: PetscInt p;
1144: if (s->maxDof == PETSC_MIN_INT) {
1145: s->maxDof = 0;
1146: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1147: s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1148: }
1149: }
1150: *maxDof = s->maxDof;
1151: return 0;
1152: }
1154: /*@
1155: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
1157: Not Collective
1159: Input Parameter:
1160: . s - the PetscSection
1162: Output Parameter:
1163: . size - the size of an array which can hold all the dofs
1165: Level: intermediate
1167: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1168: @*/
1169: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1170: {
1171: PetscInt p, n = 0;
1175: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1176: *size = n;
1177: return 0;
1178: }
1180: /*@
1181: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1183: Not Collective
1185: Input Parameter:
1186: . s - the PetscSection
1188: Output Parameter:
1189: . size - the size of an array which can hold all unconstrained dofs
1191: Level: intermediate
1193: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1194: @*/
1195: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1196: {
1197: PetscInt p, n = 0;
1201: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1202: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1203: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1204: }
1205: *size = n;
1206: return 0;
1207: }
1209: /*@
1210: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1211: the local section and an SF describing the section point overlap.
1213: Input Parameters:
1214: + s - The PetscSection for the local field layout
1215: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1216: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1217: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1219: Output Parameter:
1220: . gsection - The PetscSection for the global field layout
1222: Note: This gives negative sizes and offsets to points not owned by this process
1224: Level: intermediate
1226: .seealso: PetscSectionCreate()
1227: @*/
1228: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1229: {
1230: PetscSection gs;
1231: const PetscInt *pind = NULL;
1232: PetscInt *recv = NULL, *neg = NULL;
1233: PetscInt pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1234: PetscInt numFields, f, numComponents;
1242: PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1243: PetscSectionGetNumFields(s, &numFields);
1244: if (numFields > 0) PetscSectionSetNumFields(gs, numFields);
1245: PetscSectionGetChart(s, &pStart, &pEnd);
1246: PetscSectionSetChart(gs, pStart, pEnd);
1247: gs->includesConstraints = includeConstraints;
1248: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1249: nlocal = nroots; /* The local/leaf space matches global/root space */
1250: /* Must allocate for all points visible to SF, which may be more than this section */
1251: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1252: PetscSFGetLeafRange(sf, NULL, &maxleaf);
1255: PetscMalloc2(nroots,&neg,nlocal,&recv);
1256: PetscArrayzero(neg,nroots);
1257: }
1258: /* Mark all local points with negative dof */
1259: for (p = pStart; p < pEnd; ++p) {
1260: PetscSectionGetDof(s, p, &dof);
1261: PetscSectionSetDof(gs, p, dof);
1262: PetscSectionGetConstraintDof(s, p, &cdof);
1263: if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(gs, p, cdof);
1264: if (neg) neg[p] = -(dof+1);
1265: }
1266: PetscSectionSetUpBC(gs);
1267: if (gs->bcIndices) PetscArraycpy(gs->bcIndices, s->bcIndices, gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]);
1268: if (nroots >= 0) {
1269: PetscArrayzero(recv,nlocal);
1270: PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1271: PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1272: for (p = pStart; p < pEnd; ++p) {
1273: if (recv[p] < 0) {
1274: gs->atlasDof[p-pStart] = recv[p];
1275: PetscSectionGetDof(s, p, &dof);
1277: }
1278: }
1279: }
1280: /* Calculate new sizes, get process offset, and calculate point offsets */
1281: if (s->perm) ISGetIndices(s->perm, &pind);
1282: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1283: const PetscInt q = pind ? pind[p] : p;
1285: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1286: gs->atlasOff[q] = off;
1287: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1288: }
1289: if (!localOffsets) {
1290: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1291: globalOff -= off;
1292: }
1293: for (p = pStart, off = 0; p < pEnd; ++p) {
1294: gs->atlasOff[p-pStart] += globalOff;
1295: if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1296: }
1297: if (s->perm) ISRestoreIndices(s->perm, &pind);
1298: /* Put in negative offsets for ghost points */
1299: if (nroots >= 0) {
1300: PetscArrayzero(recv,nlocal);
1301: PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1302: PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1303: for (p = pStart; p < pEnd; ++p) {
1304: if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1305: }
1306: }
1307: PetscFree2(neg,recv);
1308: /* Set field dofs/offsets/constraints */
1309: for (f = 0; f < numFields; ++f) {
1310: gs->field[f]->includesConstraints = includeConstraints;
1311: PetscSectionGetFieldComponents(s, f, &numComponents);
1312: PetscSectionSetFieldComponents(gs, f, numComponents);
1313: }
1314: for (p = pStart; p < pEnd; ++p) {
1315: PetscSectionGetOffset(gs, p, &off);
1316: for (f = 0, foff = off; f < numFields; ++f) {
1317: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1318: if (!includeConstraints && cdof > 0) PetscSectionSetFieldConstraintDof(gs, p, f, cdof);
1319: PetscSectionGetFieldDof(s, p, f, &dof);
1320: PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof);
1321: PetscSectionSetFieldOffset(gs, p, f, foff);
1322: PetscSectionGetFieldConstraintDof(gs, p, f, &cdof);
1323: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1324: }
1325: }
1326: for (f = 0; f < numFields; ++f) {
1327: PetscSection gfs = gs->field[f];
1329: PetscSectionSetUpBC(gfs);
1330: if (gfs->bcIndices) PetscArraycpy(gfs->bcIndices, s->field[f]->bcIndices, gfs->bc->atlasOff[gfs->bc->pEnd-gfs->bc->pStart-1] + gfs->bc->atlasDof[gfs->bc->pEnd-gfs->bc->pStart-1]);
1331: }
1332: gs->setup = PETSC_TRUE;
1333: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1334: *gsection = gs;
1335: return 0;
1336: }
1338: /*@
1339: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1340: the local section and an SF describing the section point overlap.
1342: Input Parameters:
1343: + s - The PetscSection for the local field layout
1344: . sf - The SF describing parallel layout of the section points
1345: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1346: . numExcludes - The number of exclusion ranges
1347: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1349: Output Parameter:
1350: . gsection - The PetscSection for the global field layout
1352: Note: This gives negative sizes and offsets to points not owned by this process
1354: Level: advanced
1356: .seealso: PetscSectionCreate()
1357: @*/
1358: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1359: {
1360: const PetscInt *pind = NULL;
1361: PetscInt *neg = NULL, *tmpOff = NULL;
1362: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1367: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1368: PetscSectionGetChart(s, &pStart, &pEnd);
1369: PetscSectionSetChart(*gsection, pStart, pEnd);
1370: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1371: if (nroots >= 0) {
1373: PetscCalloc1(nroots, &neg);
1374: if (nroots > pEnd-pStart) {
1375: PetscCalloc1(nroots, &tmpOff);
1376: } else {
1377: tmpOff = &(*gsection)->atlasDof[-pStart];
1378: }
1379: }
1380: /* Mark ghost points with negative dof */
1381: for (p = pStart; p < pEnd; ++p) {
1382: for (e = 0; e < numExcludes; ++e) {
1383: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1384: PetscSectionSetDof(*gsection, p, 0);
1385: break;
1386: }
1387: }
1388: if (e < numExcludes) continue;
1389: PetscSectionGetDof(s, p, &dof);
1390: PetscSectionSetDof(*gsection, p, dof);
1391: PetscSectionGetConstraintDof(s, p, &cdof);
1392: if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(*gsection, p, cdof);
1393: if (neg) neg[p] = -(dof+1);
1394: }
1395: PetscSectionSetUpBC(*gsection);
1396: if (nroots >= 0) {
1397: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1398: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1399: if (nroots > pEnd - pStart) {
1400: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1401: }
1402: }
1403: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1404: if (s->perm) ISGetIndices(s->perm, &pind);
1405: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1406: const PetscInt q = pind ? pind[p] : p;
1408: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1409: (*gsection)->atlasOff[q] = off;
1410: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1411: }
1412: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1413: globalOff -= off;
1414: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1415: (*gsection)->atlasOff[p] += globalOff;
1416: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1417: }
1418: if (s->perm) ISRestoreIndices(s->perm, &pind);
1419: /* Put in negative offsets for ghost points */
1420: if (nroots >= 0) {
1421: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1422: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1423: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1424: if (nroots > pEnd - pStart) {
1425: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1426: }
1427: }
1428: if (nroots >= 0 && nroots > pEnd-pStart) PetscFree(tmpOff);
1429: PetscFree(neg);
1430: return 0;
1431: }
1433: /*@
1434: PetscSectionGetPointLayout - Get the PetscLayout associated with the section points
1436: Collective on comm
1438: Input Parameters:
1439: + comm - The MPI_Comm
1440: - s - The PetscSection
1442: Output Parameter:
1443: . layout - The point layout for the section
1445: Note: This is usually called for the default global section.
1447: Level: advanced
1449: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1450: @*/
1451: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1452: {
1453: PetscInt pStart, pEnd, p, localSize = 0;
1455: PetscSectionGetChart(s, &pStart, &pEnd);
1456: for (p = pStart; p < pEnd; ++p) {
1457: PetscInt dof;
1459: PetscSectionGetDof(s, p, &dof);
1460: if (dof >= 0) ++localSize;
1461: }
1462: PetscLayoutCreate(comm, layout);
1463: PetscLayoutSetLocalSize(*layout, localSize);
1464: PetscLayoutSetBlockSize(*layout, 1);
1465: PetscLayoutSetUp(*layout);
1466: return 0;
1467: }
1469: /*@
1470: PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.
1472: Collective on comm
1474: Input Parameters:
1475: + comm - The MPI_Comm
1476: - s - The PetscSection
1478: Output Parameter:
1479: . layout - The dof layout for the section
1481: Note: This is usually called for the default global section.
1483: Level: advanced
1485: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1486: @*/
1487: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1488: {
1489: PetscInt pStart, pEnd, p, localSize = 0;
1493: PetscSectionGetChart(s, &pStart, &pEnd);
1494: for (p = pStart; p < pEnd; ++p) {
1495: PetscInt dof,cdof;
1497: PetscSectionGetDof(s, p, &dof);
1498: PetscSectionGetConstraintDof(s, p, &cdof);
1499: if (dof-cdof > 0) localSize += dof-cdof;
1500: }
1501: PetscLayoutCreate(comm, layout);
1502: PetscLayoutSetLocalSize(*layout, localSize);
1503: PetscLayoutSetBlockSize(*layout, 1);
1504: PetscLayoutSetUp(*layout);
1505: return 0;
1506: }
1508: /*@
1509: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1511: Not Collective
1513: Input Parameters:
1514: + s - the PetscSection
1515: - point - the point
1517: Output Parameter:
1518: . offset - the offset
1520: Level: intermediate
1522: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1523: @*/
1524: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1525: {
1528: if (PetscDefined(USE_DEBUG)) {
1530: }
1531: *offset = s->atlasOff[point - s->pStart];
1532: return 0;
1533: }
1535: /*@
1536: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1538: Not Collective
1540: Input Parameters:
1541: + s - the PetscSection
1542: . point - the point
1543: - offset - the offset
1545: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1547: Level: intermediate
1549: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1550: @*/
1551: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1552: {
1555: s->atlasOff[point - s->pStart] = offset;
1556: return 0;
1557: }
1559: /*@
1560: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1562: Not Collective
1564: Input Parameters:
1565: + s - the PetscSection
1566: . point - the point
1567: - field - the field
1569: Output Parameter:
1570: . offset - the offset
1572: Level: intermediate
1574: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1575: @*/
1576: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1577: {
1580: PetscSectionCheckValidField(field,s->numFields);
1581: PetscSectionGetOffset(s->field[field], point, offset);
1582: return 0;
1583: }
1585: /*@
1586: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given field at a point.
1588: Not Collective
1590: Input Parameters:
1591: + s - the PetscSection
1592: . point - the point
1593: . field - the field
1594: - offset - the offset
1596: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1598: Level: intermediate
1600: .seealso: PetscSectionGetFieldOffset(), PetscSectionSetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1601: @*/
1602: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1603: {
1605: PetscSectionCheckValidField(field,s->numFields);
1606: PetscSectionSetOffset(s->field[field], point, offset);
1607: return 0;
1608: }
1610: /*@
1611: PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.
1613: Not Collective
1615: Input Parameters:
1616: + s - the PetscSection
1617: . point - the point
1618: - field - the field
1620: Output Parameter:
1621: . offset - the offset
1623: Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1624: this point, what number is the first dof with this field.
1626: Level: advanced
1628: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1629: @*/
1630: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1631: {
1632: PetscInt off, foff;
1636: PetscSectionCheckValidField(field,s->numFields);
1637: PetscSectionGetOffset(s, point, &off);
1638: PetscSectionGetOffset(s->field[field], point, &foff);
1639: *offset = foff - off;
1640: return 0;
1641: }
1643: /*@
1644: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1646: Not Collective
1648: Input Parameter:
1649: . s - the PetscSection
1651: Output Parameters:
1652: + start - the minimum offset
1653: - end - one more than the maximum offset
1655: Level: intermediate
1657: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1658: @*/
1659: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1660: {
1661: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1664: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1665: PetscSectionGetChart(s, &pStart, &pEnd);
1666: for (p = 0; p < pEnd-pStart; ++p) {
1667: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1669: if (off >= 0) {
1670: os = PetscMin(os, off);
1671: oe = PetscMax(oe, off+dof);
1672: }
1673: }
1674: if (start) *start = os;
1675: if (end) *end = oe;
1676: return 0;
1677: }
1679: /*@
1680: PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields
1682: Collective
1684: Input Parameters:
1685: + s - the PetscSection
1686: . len - the number of subfields
1687: - fields - the subfield numbers
1689: Output Parameter:
1690: . subs - the subsection
1692: Note: The section offsets now refer to a new, smaller vector.
1694: Level: advanced
1696: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1697: @*/
1698: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1699: {
1700: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
1702: if (!len) return 0;
1706: PetscSectionGetNumFields(s, &nF);
1708: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1709: PetscSectionSetNumFields(*subs, len);
1710: for (f = 0; f < len; ++f) {
1711: const char *name = NULL;
1712: PetscInt numComp = 0;
1714: PetscSectionGetFieldName(s, fields[f], &name);
1715: PetscSectionSetFieldName(*subs, f, name);
1716: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1717: PetscSectionSetFieldComponents(*subs, f, numComp);
1718: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1719: PetscSectionGetComponentName(s, fields[f], c, &name);
1720: PetscSectionSetComponentName(*subs, f, c, name);
1721: }
1722: }
1723: PetscSectionGetChart(s, &pStart, &pEnd);
1724: PetscSectionSetChart(*subs, pStart, pEnd);
1725: for (p = pStart; p < pEnd; ++p) {
1726: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1728: for (f = 0; f < len; ++f) {
1729: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1730: PetscSectionSetFieldDof(*subs, p, f, fdof);
1731: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1732: if (cfdof) PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);
1733: dof += fdof;
1734: cdof += cfdof;
1735: }
1736: PetscSectionSetDof(*subs, p, dof);
1737: if (cdof) PetscSectionSetConstraintDof(*subs, p, cdof);
1738: maxCdof = PetscMax(cdof, maxCdof);
1739: }
1740: PetscSectionSetUp(*subs);
1741: if (maxCdof) {
1742: PetscInt *indices;
1744: PetscMalloc1(maxCdof, &indices);
1745: for (p = pStart; p < pEnd; ++p) {
1746: PetscInt cdof;
1748: PetscSectionGetConstraintDof(*subs, p, &cdof);
1749: if (cdof) {
1750: const PetscInt *oldIndices = NULL;
1751: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1753: for (f = 0; f < len; ++f) {
1754: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1755: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1756: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1757: PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices);
1758: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1759: numConst += cfdof;
1760: fOff += fdof;
1761: }
1762: PetscSectionSetConstraintIndices(*subs, p, indices);
1763: }
1764: }
1765: PetscFree(indices);
1766: }
1767: return 0;
1768: }
1770: /*@
1771: PetscSectionCreateSupersection - Create a new, larger section composed of the input sections
1773: Collective
1775: Input Parameters:
1776: + s - the input sections
1777: - len - the number of input sections
1779: Output Parameter:
1780: . supers - the supersection
1782: Note: The section offsets now refer to a new, larger vector.
1784: Level: advanced
1786: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1787: @*/
1788: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1789: {
1790: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1792: if (!len) return 0;
1793: for (i = 0; i < len; ++i) {
1794: PetscInt nf, pStarti, pEndi;
1796: PetscSectionGetNumFields(s[i], &nf);
1797: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1798: pStart = PetscMin(pStart, pStarti);
1799: pEnd = PetscMax(pEnd, pEndi);
1800: Nf += nf;
1801: }
1802: PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1803: PetscSectionSetNumFields(*supers, Nf);
1804: for (i = 0, f = 0; i < len; ++i) {
1805: PetscInt nf, fi, ci;
1807: PetscSectionGetNumFields(s[i], &nf);
1808: for (fi = 0; fi < nf; ++fi, ++f) {
1809: const char *name = NULL;
1810: PetscInt numComp = 0;
1812: PetscSectionGetFieldName(s[i], fi, &name);
1813: PetscSectionSetFieldName(*supers, f, name);
1814: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1815: PetscSectionSetFieldComponents(*supers, f, numComp);
1816: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1817: PetscSectionGetComponentName(s[i], fi, ci, &name);
1818: PetscSectionSetComponentName(*supers, f, ci, name);
1819: }
1820: }
1821: }
1822: PetscSectionSetChart(*supers, pStart, pEnd);
1823: for (p = pStart; p < pEnd; ++p) {
1824: PetscInt dof = 0, cdof = 0;
1826: for (i = 0, f = 0; i < len; ++i) {
1827: PetscInt nf, fi, pStarti, pEndi;
1828: PetscInt fdof = 0, cfdof = 0;
1830: PetscSectionGetNumFields(s[i], &nf);
1831: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1832: if ((p < pStarti) || (p >= pEndi)) continue;
1833: for (fi = 0; fi < nf; ++fi, ++f) {
1834: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1835: PetscSectionAddFieldDof(*supers, p, f, fdof);
1836: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1837: if (cfdof) PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);
1838: dof += fdof;
1839: cdof += cfdof;
1840: }
1841: }
1842: PetscSectionSetDof(*supers, p, dof);
1843: if (cdof) PetscSectionSetConstraintDof(*supers, p, cdof);
1844: maxCdof = PetscMax(cdof, maxCdof);
1845: }
1846: PetscSectionSetUp(*supers);
1847: if (maxCdof) {
1848: PetscInt *indices;
1850: PetscMalloc1(maxCdof, &indices);
1851: for (p = pStart; p < pEnd; ++p) {
1852: PetscInt cdof;
1854: PetscSectionGetConstraintDof(*supers, p, &cdof);
1855: if (cdof) {
1856: PetscInt dof, numConst = 0, fOff = 0;
1858: for (i = 0, f = 0; i < len; ++i) {
1859: const PetscInt *oldIndices = NULL;
1860: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1862: PetscSectionGetNumFields(s[i], &nf);
1863: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1864: if ((p < pStarti) || (p >= pEndi)) continue;
1865: for (fi = 0; fi < nf; ++fi, ++f) {
1866: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1867: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1868: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1869: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc];
1870: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1871: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] += fOff;
1872: numConst += cfdof;
1873: }
1874: PetscSectionGetDof(s[i], p, &dof);
1875: fOff += dof;
1876: }
1877: PetscSectionSetConstraintIndices(*supers, p, indices);
1878: }
1879: }
1880: PetscFree(indices);
1881: }
1882: return 0;
1883: }
1885: /*@
1886: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
1888: Collective
1890: Input Parameters:
1891: + s - the PetscSection
1892: - subpointMap - a sorted list of points in the original mesh which are in the submesh
1894: Output Parameter:
1895: . subs - the subsection
1897: Note: The section offsets now refer to a new, smaller vector.
1899: Level: advanced
1901: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1902: @*/
1903: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1904: {
1905: const PetscInt *points = NULL, *indices = NULL;
1906: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;
1911: PetscSectionGetNumFields(s, &numFields);
1912: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1913: if (numFields) PetscSectionSetNumFields(*subs, numFields);
1914: for (f = 0; f < numFields; ++f) {
1915: const char *name = NULL;
1916: PetscInt numComp = 0;
1918: PetscSectionGetFieldName(s, f, &name);
1919: PetscSectionSetFieldName(*subs, f, name);
1920: PetscSectionGetFieldComponents(s, f, &numComp);
1921: PetscSectionSetFieldComponents(*subs, f, numComp);
1922: for (c = 0; c < s->numFieldComponents[f]; ++c) {
1923: PetscSectionGetComponentName(s, f, c, &name);
1924: PetscSectionSetComponentName(*subs, f, c, name);
1925: }
1926: }
1927: /* For right now, we do not try to squeeze the subchart */
1928: if (subpointMap) {
1929: ISGetSize(subpointMap, &numSubpoints);
1930: ISGetIndices(subpointMap, &points);
1931: }
1932: PetscSectionGetChart(s, &pStart, &pEnd);
1933: PetscSectionSetChart(*subs, 0, numSubpoints);
1934: for (p = pStart; p < pEnd; ++p) {
1935: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1937: PetscFindInt(p, numSubpoints, points, &subp);
1938: if (subp < 0) continue;
1939: for (f = 0; f < numFields; ++f) {
1940: PetscSectionGetFieldDof(s, p, f, &fdof);
1941: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1942: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1943: if (cfdof) PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);
1944: }
1945: PetscSectionGetDof(s, p, &dof);
1946: PetscSectionSetDof(*subs, subp, dof);
1947: PetscSectionGetConstraintDof(s, p, &cdof);
1948: if (cdof) PetscSectionSetConstraintDof(*subs, subp, cdof);
1949: }
1950: PetscSectionSetUp(*subs);
1951: /* Change offsets to original offsets */
1952: for (p = pStart; p < pEnd; ++p) {
1953: PetscInt off, foff = 0;
1955: PetscFindInt(p, numSubpoints, points, &subp);
1956: if (subp < 0) continue;
1957: for (f = 0; f < numFields; ++f) {
1958: PetscSectionGetFieldOffset(s, p, f, &foff);
1959: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1960: }
1961: PetscSectionGetOffset(s, p, &off);
1962: PetscSectionSetOffset(*subs, subp, off);
1963: }
1964: /* Copy constraint indices */
1965: for (subp = 0; subp < numSubpoints; ++subp) {
1966: PetscInt cdof;
1968: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1969: if (cdof) {
1970: for (f = 0; f < numFields; ++f) {
1971: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1972: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1973: }
1974: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1975: PetscSectionSetConstraintIndices(*subs, subp, indices);
1976: }
1977: }
1978: if (subpointMap) ISRestoreIndices(subpointMap, &points);
1979: return 0;
1980: }
1982: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1983: {
1984: PetscInt p;
1985: PetscMPIInt rank;
1987: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1988: PetscViewerASCIIPushSynchronized(viewer);
1989: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1990: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1991: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1992: PetscInt b;
1994: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1995: if (s->bcIndices) {
1996: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1997: PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p]+b]);
1998: }
1999: }
2000: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2001: } else {
2002: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2003: }
2004: }
2005: PetscViewerFlush(viewer);
2006: PetscViewerASCIIPopSynchronized(viewer);
2007: if (s->sym) {
2008: PetscViewerASCIIPushTab(viewer);
2009: PetscSectionSymView(s->sym,viewer);
2010: PetscViewerASCIIPopTab(viewer);
2011: }
2012: return 0;
2013: }
2015: /*@C
2016: PetscSectionViewFromOptions - View from Options
2018: Collective
2020: Input Parameters:
2021: + A - the PetscSection object to view
2022: . obj - Optional object
2023: - name - command line option
2025: Level: intermediate
2026: .seealso: PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
2027: @*/
2028: PetscErrorCode PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
2029: {
2031: PetscObjectViewFromOptions((PetscObject)A,obj,name);
2032: return 0;
2033: }
2035: /*@C
2036: PetscSectionView - Views a PetscSection
2038: Collective
2040: Input Parameters:
2041: + s - the PetscSection object to view
2042: - v - the viewer
2044: Note:
2045: PetscSectionView(), when viewer is of type PETSCVIEWERHDF5, only saves
2046: distribution independent data, such as dofs, offsets, constraint dofs,
2047: and constraint indices. Points that have negative dofs, for instance,
2048: are not saved as they represent points owned by other processes.
2049: Point numbering and rank assignment is currently not stored.
2050: The saved section can be loaded with PetscSectionLoad().
2052: Level: beginner
2054: .seealso PetscSectionCreate(), PetscSectionDestroy(), PetscSectionLoad()
2055: @*/
2056: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2057: {
2058: PetscBool isascii, ishdf5;
2059: PetscInt f;
2062: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);
2064: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
2065: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
2066: if (isascii) {
2067: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
2068: if (s->numFields) {
2069: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields);
2070: for (f = 0; f < s->numFields; ++f) {
2071: PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
2072: PetscSectionView_ASCII(s->field[f], viewer);
2073: }
2074: } else {
2075: PetscSectionView_ASCII(s, viewer);
2076: }
2077: } else if (ishdf5) {
2078: #if PetscDefined(HAVE_HDF5)
2079: PetscSectionView_HDF5_Internal(s, viewer);
2080: #else
2081: SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2082: #endif
2083: }
2084: return 0;
2085: }
2087: /*@C
2088: PetscSectionLoad - Loads a PetscSection
2090: Collective
2092: Input Parameters:
2093: + s - the PetscSection object to load
2094: - v - the viewer
2096: Note:
2097: PetscSectionLoad(), when viewer is of type PETSCVIEWERHDF5, loads
2098: a section saved with PetscSectionView(). The number of processes
2099: used here (N) does not need to be the same as that used when saving.
2100: After calling this function, the chart of s on rank i will be set
2101: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2102: saved section points.
2104: Level: beginner
2106: .seealso PetscSectionCreate(), PetscSectionDestroy(), PetscSectionView()
2107: @*/
2108: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2109: {
2110: PetscBool ishdf5;
2114: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
2115: if (ishdf5) {
2116: #if PetscDefined(HAVE_HDF5)
2117: PetscSectionLoad_HDF5_Internal(s, viewer);
2118: return 0;
2119: #else
2120: SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2121: #endif
2122: } else SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2123: }
2125: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2126: {
2127: PetscSectionClosurePermVal clVal;
2129: if (!section->clHash) return 0;
2130: kh_foreach_value(section->clHash, clVal, {
2131: PetscFree(clVal.perm);
2132: PetscFree(clVal.invPerm);
2133: });
2134: kh_destroy(ClPerm, section->clHash);
2135: section->clHash = NULL;
2136: return 0;
2137: }
2139: /*@
2140: PetscSectionReset - Frees all section data.
2142: Not Collective
2144: Input Parameters:
2145: . s - the PetscSection
2147: Level: beginner
2149: .seealso: PetscSection, PetscSectionCreate()
2150: @*/
2151: PetscErrorCode PetscSectionReset(PetscSection s)
2152: {
2153: PetscInt f, c;
2156: for (f = 0; f < s->numFields; ++f) {
2157: PetscSectionDestroy(&s->field[f]);
2158: PetscFree(s->fieldNames[f]);
2159: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2160: PetscFree(s->compNames[f][c]);
2161: }
2162: PetscFree(s->compNames[f]);
2163: }
2164: PetscFree(s->numFieldComponents);
2165: PetscFree(s->fieldNames);
2166: PetscFree(s->compNames);
2167: PetscFree(s->field);
2168: PetscSectionDestroy(&s->bc);
2169: PetscFree(s->bcIndices);
2170: PetscFree2(s->atlasDof, s->atlasOff);
2171: PetscSectionDestroy(&s->clSection);
2172: ISDestroy(&s->clPoints);
2173: ISDestroy(&s->perm);
2174: PetscSectionResetClosurePermutation(s);
2175: PetscSectionSymDestroy(&s->sym);
2176: PetscSectionDestroy(&s->clSection);
2177: ISDestroy(&s->clPoints);
2178: PetscSectionInvalidateMaxDof_Internal(s);
2179: s->pStart = -1;
2180: s->pEnd = -1;
2181: s->maxDof = 0;
2182: s->setup = PETSC_FALSE;
2183: s->numFields = 0;
2184: s->clObj = NULL;
2185: return 0;
2186: }
2188: /*@
2189: PetscSectionDestroy - Frees a section object and frees its range if that exists.
2191: Not Collective
2193: Input Parameters:
2194: . s - the PetscSection
2196: Level: beginner
2198: .seealso: PetscSection, PetscSectionCreate()
2199: @*/
2200: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2201: {
2202: if (!*s) return 0;
2204: if (--((PetscObject)(*s))->refct > 0) {
2205: *s = NULL;
2206: return 0;
2207: }
2208: PetscSectionReset(*s);
2209: PetscHeaderDestroy(s);
2210: return 0;
2211: }
2213: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2214: {
2215: const PetscInt p = point - s->pStart;
2218: *values = &baseArray[s->atlasOff[p]];
2219: return 0;
2220: }
2222: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2223: {
2224: PetscInt *array;
2225: const PetscInt p = point - s->pStart;
2226: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2227: PetscInt cDim = 0;
2230: PetscSectionGetConstraintDof(s, p, &cDim);
2231: array = &baseArray[s->atlasOff[p]];
2232: if (!cDim) {
2233: if (orientation >= 0) {
2234: const PetscInt dim = s->atlasDof[p];
2235: PetscInt i;
2237: if (mode == INSERT_VALUES) {
2238: for (i = 0; i < dim; ++i) array[i] = values[i];
2239: } else {
2240: for (i = 0; i < dim; ++i) array[i] += values[i];
2241: }
2242: } else {
2243: PetscInt offset = 0;
2244: PetscInt j = -1, field, i;
2246: for (field = 0; field < s->numFields; ++field) {
2247: const PetscInt dim = s->field[field]->atlasDof[p];
2249: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2250: offset += dim;
2251: }
2252: }
2253: } else {
2254: if (orientation >= 0) {
2255: const PetscInt dim = s->atlasDof[p];
2256: PetscInt cInd = 0, i;
2257: const PetscInt *cDof;
2259: PetscSectionGetConstraintIndices(s, point, &cDof);
2260: if (mode == INSERT_VALUES) {
2261: for (i = 0; i < dim; ++i) {
2262: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2263: array[i] = values[i];
2264: }
2265: } else {
2266: for (i = 0; i < dim; ++i) {
2267: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2268: array[i] += values[i];
2269: }
2270: }
2271: } else {
2272: const PetscInt *cDof;
2273: PetscInt offset = 0;
2274: PetscInt cOffset = 0;
2275: PetscInt j = 0, field;
2277: PetscSectionGetConstraintIndices(s, point, &cDof);
2278: for (field = 0; field < s->numFields; ++field) {
2279: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2280: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2281: const PetscInt sDim = dim - tDim;
2282: PetscInt cInd = 0, i,k;
2284: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2285: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2286: array[j] = values[k];
2287: }
2288: offset += dim;
2289: cOffset += dim - tDim;
2290: }
2291: }
2292: }
2293: return 0;
2294: }
2296: /*@C
2297: PetscSectionHasConstraints - Determine whether a section has constrained dofs
2299: Not Collective
2301: Input Parameter:
2302: . s - The PetscSection
2304: Output Parameter:
2305: . hasConstraints - flag indicating that the section has constrained dofs
2307: Level: intermediate
2309: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2310: @*/
2311: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2312: {
2315: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2316: return 0;
2317: }
2319: /*@C
2320: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2322: Not Collective
2324: Input Parameters:
2325: + s - The PetscSection
2326: - point - The point
2328: Output Parameter:
2329: . indices - The constrained dofs
2331: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
2333: Level: intermediate
2335: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2336: @*/
2337: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2338: {
2340: if (s->bc) {
2341: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2342: } else *indices = NULL;
2343: return 0;
2344: }
2346: /*@C
2347: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2349: Not Collective
2351: Input Parameters:
2352: + s - The PetscSection
2353: . point - The point
2354: - indices - The constrained dofs
2356: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
2358: Level: intermediate
2360: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2361: @*/
2362: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2363: {
2365: if (s->bc) {
2366: const PetscInt dof = s->atlasDof[point];
2367: const PetscInt cdof = s->bc->atlasDof[point];
2368: PetscInt d;
2370: for (d = 0; d < cdof; ++d) {
2372: }
2373: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2374: }
2375: return 0;
2376: }
2378: /*@C
2379: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2381: Not Collective
2383: Input Parameters:
2384: + s - The PetscSection
2385: . field - The field number
2386: - point - The point
2388: Output Parameter:
2389: . indices - The constrained dofs sorted in ascending order
2391: Notes:
2392: The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by PetscSectionGetConstraintDof().
2394: Fortran Note:
2395: In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()
2397: Level: intermediate
2399: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2400: @*/
2401: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2402: {
2405: PetscSectionCheckValidField(field,s->numFields);
2406: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2407: return 0;
2408: }
2410: /*@C
2411: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2413: Not Collective
2415: Input Parameters:
2416: + s - The PetscSection
2417: . point - The point
2418: . field - The field number
2419: - indices - The constrained dofs
2421: Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()
2423: Level: intermediate
2425: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2426: @*/
2427: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2428: {
2430: if (PetscDefined(USE_DEBUG)) {
2431: PetscInt nfdof;
2433: PetscSectionGetFieldConstraintDof(s, point, field, &nfdof);
2435: }
2436: PetscSectionCheckValidField(field,s->numFields);
2437: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2438: return 0;
2439: }
2441: /*@
2442: PetscSectionPermute - Reorder the section according to the input point permutation
2444: Collective
2446: Input Parameters:
2447: + section - The PetscSection object
2448: - perm - The point permutation, old point p becomes new point perm[p]
2450: Output Parameter:
2451: . sectionNew - The permuted PetscSection
2453: Level: intermediate
2455: .seealso: MatPermute()
2456: @*/
2457: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2458: {
2459: PetscSection s = section, sNew;
2460: const PetscInt *perm;
2461: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2466: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2467: PetscSectionGetNumFields(s, &numFields);
2468: if (numFields) PetscSectionSetNumFields(sNew, numFields);
2469: for (f = 0; f < numFields; ++f) {
2470: const char *name;
2471: PetscInt numComp;
2473: PetscSectionGetFieldName(s, f, &name);
2474: PetscSectionSetFieldName(sNew, f, name);
2475: PetscSectionGetFieldComponents(s, f, &numComp);
2476: PetscSectionSetFieldComponents(sNew, f, numComp);
2477: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2478: PetscSectionGetComponentName(s, f, c, &name);
2479: PetscSectionSetComponentName(sNew, f, c, name);
2480: }
2481: }
2482: ISGetLocalSize(permutation, &numPoints);
2483: ISGetIndices(permutation, &perm);
2484: PetscSectionGetChart(s, &pStart, &pEnd);
2485: PetscSectionSetChart(sNew, pStart, pEnd);
2487: for (p = pStart; p < pEnd; ++p) {
2488: PetscInt dof, cdof;
2490: PetscSectionGetDof(s, p, &dof);
2491: PetscSectionSetDof(sNew, perm[p], dof);
2492: PetscSectionGetConstraintDof(s, p, &cdof);
2493: if (cdof) PetscSectionSetConstraintDof(sNew, perm[p], cdof);
2494: for (f = 0; f < numFields; ++f) {
2495: PetscSectionGetFieldDof(s, p, f, &dof);
2496: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2497: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2498: if (cdof) PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);
2499: }
2500: }
2501: PetscSectionSetUp(sNew);
2502: for (p = pStart; p < pEnd; ++p) {
2503: const PetscInt *cind;
2504: PetscInt cdof;
2506: PetscSectionGetConstraintDof(s, p, &cdof);
2507: if (cdof) {
2508: PetscSectionGetConstraintIndices(s, p, &cind);
2509: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2510: }
2511: for (f = 0; f < numFields; ++f) {
2512: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2513: if (cdof) {
2514: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2515: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2516: }
2517: }
2518: }
2519: ISRestoreIndices(permutation, &perm);
2520: *sectionNew = sNew;
2521: return 0;
2522: }
2524: /*@
2525: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2527: Collective
2529: Input Parameters:
2530: + section - The PetscSection
2531: . obj - A PetscObject which serves as the key for this index
2532: . clSection - Section giving the size of the closure of each point
2533: - clPoints - IS giving the points in each closure
2535: Note: We compress out closure points with no dofs in this section
2537: Level: advanced
2539: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2540: @*/
2541: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2542: {
2546: if (section->clObj != obj) PetscSectionResetClosurePermutation(section);
2547: section->clObj = obj;
2548: PetscObjectReference((PetscObject)clSection);
2549: PetscObjectReference((PetscObject)clPoints);
2550: PetscSectionDestroy(§ion->clSection);
2551: ISDestroy(§ion->clPoints);
2552: section->clSection = clSection;
2553: section->clPoints = clPoints;
2554: return 0;
2555: }
2557: /*@
2558: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2560: Collective
2562: Input Parameters:
2563: + section - The PetscSection
2564: - obj - A PetscObject which serves as the key for this index
2566: Output Parameters:
2567: + clSection - Section giving the size of the closure of each point
2568: - clPoints - IS giving the points in each closure
2570: Note: We compress out closure points with no dofs in this section
2572: Level: advanced
2574: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2575: @*/
2576: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2577: {
2578: if (section->clObj == obj) {
2579: if (clSection) *clSection = section->clSection;
2580: if (clPoints) *clPoints = section->clPoints;
2581: } else {
2582: if (clSection) *clSection = NULL;
2583: if (clPoints) *clPoints = NULL;
2584: }
2585: return 0;
2586: }
2588: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2589: {
2590: PetscInt i;
2591: khiter_t iter;
2592: int new_entry;
2593: PetscSectionClosurePermKey key = {depth, clSize};
2594: PetscSectionClosurePermVal *val;
2596: if (section->clObj != obj) {
2597: PetscSectionDestroy(§ion->clSection);
2598: ISDestroy(§ion->clPoints);
2599: }
2600: section->clObj = obj;
2601: if (!section->clHash) PetscClPermCreate(§ion->clHash);
2602: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2603: val = &kh_val(section->clHash, iter);
2604: if (!new_entry) {
2605: PetscFree(val->perm);
2606: PetscFree(val->invPerm);
2607: }
2608: if (mode == PETSC_COPY_VALUES) {
2609: PetscMalloc1(clSize, &val->perm);
2610: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2611: PetscArraycpy(val->perm, clPerm, clSize);
2612: } else if (mode == PETSC_OWN_POINTER) {
2613: val->perm = clPerm;
2614: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2615: PetscMalloc1(clSize, &val->invPerm);
2616: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2617: return 0;
2618: }
2620: /*@
2621: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2623: Not Collective
2625: Input Parameters:
2626: + section - The PetscSection
2627: . obj - A PetscObject which serves as the key for this index (usually a DM)
2628: . depth - Depth of points on which to apply the given permutation
2629: - perm - Permutation of the cell dof closure
2631: Note:
2632: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2633: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2634: each topology and degree.
2636: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2637: exotic/enriched spaces on mixed topology meshes.
2639: Level: intermediate
2641: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2642: @*/
2643: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2644: {
2645: const PetscInt *clPerm = NULL;
2646: PetscInt clSize = 0;
2648: if (perm) {
2649: ISGetLocalSize(perm, &clSize);
2650: ISGetIndices(perm, &clPerm);
2651: }
2652: PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2653: if (perm) ISRestoreIndices(perm, &clPerm);
2654: return 0;
2655: }
2657: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2658: {
2659: if (section->clObj == obj) {
2660: PetscSectionClosurePermKey k = {depth, size};
2661: PetscSectionClosurePermVal v;
2662: PetscClPermGet(section->clHash, k, &v);
2663: if (perm) *perm = v.perm;
2664: } else {
2665: if (perm) *perm = NULL;
2666: }
2667: return 0;
2668: }
2670: /*@
2671: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2673: Not Collective
2675: Input Parameters:
2676: + section - The PetscSection
2677: . obj - A PetscObject which serves as the key for this index (usually a DM)
2678: . depth - Depth stratum on which to obtain closure permutation
2679: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2681: Output Parameter:
2682: . perm - The dof closure permutation
2684: Note:
2685: The user must destroy the IS that is returned.
2687: Level: intermediate
2689: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2690: @*/
2691: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2692: {
2693: const PetscInt *clPerm;
2695: PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2696: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2697: return 0;
2698: }
2700: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2701: {
2702: if (section->clObj == obj && section->clHash) {
2703: PetscSectionClosurePermKey k = {depth, size};
2704: PetscSectionClosurePermVal v;
2705: PetscClPermGet(section->clHash, k, &v);
2706: if (perm) *perm = v.invPerm;
2707: } else {
2708: if (perm) *perm = NULL;
2709: }
2710: return 0;
2711: }
2713: /*@
2714: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2716: Not Collective
2718: Input Parameters:
2719: + section - The PetscSection
2720: . obj - A PetscObject which serves as the key for this index (usually a DM)
2721: . depth - Depth stratum on which to obtain closure permutation
2722: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2724: Output Parameters:
2725: . perm - The dof closure permutation
2727: Note:
2728: The user must destroy the IS that is returned.
2730: Level: intermediate
2732: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2733: @*/
2734: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2735: {
2736: const PetscInt *clPerm;
2738: PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2739: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2740: return 0;
2741: }
2743: /*@
2744: PetscSectionGetField - Get the subsection associated with a single field
2746: Input Parameters:
2747: + s - The PetscSection
2748: - field - The field number
2750: Output Parameter:
2751: . subs - The subsection for the given field
2753: Level: intermediate
2755: .seealso: PetscSectionSetNumFields()
2756: @*/
2757: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2758: {
2761: PetscSectionCheckValidField(field,s->numFields);
2762: *subs = s->field[field];
2763: return 0;
2764: }
2766: PetscClassId PETSC_SECTION_SYM_CLASSID;
2767: PetscFunctionList PetscSectionSymList = NULL;
2769: /*@
2770: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
2772: Collective
2774: Input Parameter:
2775: . comm - the MPI communicator
2777: Output Parameter:
2778: . sym - pointer to the new set of symmetries
2780: Level: developer
2782: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2783: @*/
2784: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2785: {
2787: ISInitializePackage();
2788: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2789: return 0;
2790: }
2792: /*@C
2793: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
2795: Collective
2797: Input Parameters:
2798: + sym - The section symmetry object
2799: - method - The name of the section symmetry type
2801: Level: developer
2803: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2804: @*/
2805: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2806: {
2807: PetscErrorCode (*r)(PetscSectionSym);
2808: PetscBool match;
2811: PetscObjectTypeCompare((PetscObject) sym, method, &match);
2812: if (match) return 0;
2814: PetscFunctionListFind(PetscSectionSymList,method,&r);
2816: if (sym->ops->destroy) {
2817: (*sym->ops->destroy)(sym);
2818: sym->ops->destroy = NULL;
2819: }
2820: (*r)(sym);
2821: PetscObjectChangeTypeName((PetscObject)sym,method);
2822: return 0;
2823: }
2825: /*@C
2826: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
2828: Not Collective
2830: Input Parameter:
2831: . sym - The section symmetry
2833: Output Parameter:
2834: . type - The index set type name
2836: Level: developer
2838: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2839: @*/
2840: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2841: {
2844: *type = ((PetscObject)sym)->type_name;
2845: return 0;
2846: }
2848: /*@C
2849: PetscSectionSymRegister - Adds a new section symmetry implementation
2851: Not Collective
2853: Input Parameters:
2854: + name - The name of a new user-defined creation routine
2855: - create_func - The creation routine itself
2857: Notes:
2858: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
2860: Level: developer
2862: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2863: @*/
2864: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2865: {
2866: ISInitializePackage();
2867: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2868: return 0;
2869: }
2871: /*@
2872: PetscSectionSymDestroy - Destroys a section symmetry.
2874: Collective
2876: Input Parameters:
2877: . sym - the section symmetry
2879: Level: developer
2881: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2882: @*/
2883: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2884: {
2885: SymWorkLink link,next;
2887: if (!*sym) return 0;
2889: if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; return 0;}
2890: if ((*sym)->ops->destroy) {
2891: (*(*sym)->ops->destroy)(*sym);
2892: }
2894: for (link=(*sym)->workin; link; link=next) {
2895: PetscInt **perms = (PetscInt **)link->perms;
2896: PetscScalar **rots = (PetscScalar **) link->rots;
2897: PetscFree2(perms,rots);
2898: next = link->next;
2899: PetscFree(link);
2900: }
2901: (*sym)->workin = NULL;
2902: PetscHeaderDestroy(sym);
2903: return 0;
2904: }
2906: /*@C
2907: PetscSectionSymView - Displays a section symmetry
2909: Collective
2911: Input Parameters:
2912: + sym - the index set
2913: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
2915: Level: developer
2917: .seealso: PetscViewerASCIIOpen()
2918: @*/
2919: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2920: {
2922: if (!viewer) {
2923: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2924: }
2927: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2928: if (sym->ops->view) {
2929: (*sym->ops->view)(sym,viewer);
2930: }
2931: return 0;
2932: }
2934: /*@
2935: PetscSectionSetSym - Set the symmetries for the data referred to by the section
2937: Collective
2939: Input Parameters:
2940: + section - the section describing data layout
2941: - sym - the symmetry describing the affect of orientation on the access of the data
2943: Level: developer
2945: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2946: @*/
2947: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2948: {
2950: PetscSectionSymDestroy(&(section->sym));
2951: if (sym) {
2954: PetscObjectReference((PetscObject) sym);
2955: }
2956: section->sym = sym;
2957: return 0;
2958: }
2960: /*@
2961: PetscSectionGetSym - Get the symmetries for the data referred to by the section
2963: Not Collective
2965: Input Parameters:
2966: . section - the section describing data layout
2968: Output Parameters:
2969: . sym - the symmetry describing the affect of orientation on the access of the data
2971: Level: developer
2973: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2974: @*/
2975: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2976: {
2978: *sym = section->sym;
2979: return 0;
2980: }
2982: /*@
2983: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
2985: Collective
2987: Input Parameters:
2988: + section - the section describing data layout
2989: . field - the field number
2990: - sym - the symmetry describing the affect of orientation on the access of the data
2992: Level: developer
2994: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2995: @*/
2996: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2997: {
2999: PetscSectionCheckValidField(field,section->numFields);
3000: PetscSectionSetSym(section->field[field],sym);
3001: return 0;
3002: }
3004: /*@
3005: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3007: Collective
3009: Input Parameters:
3010: + section - the section describing data layout
3011: - field - the field number
3013: Output Parameters:
3014: . sym - the symmetry describing the affect of orientation on the access of the data
3016: Level: developer
3018: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3019: @*/
3020: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3021: {
3023: PetscSectionCheckValidField(field,section->numFields);
3024: *sym = section->field[field]->sym;
3025: return 0;
3026: }
3028: /*@C
3029: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
3031: Not Collective
3033: Input Parameters:
3034: + section - the section
3035: . numPoints - the number of points
3036: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3037: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3038: context, see DMPlexGetConeOrientation()).
3040: Output Parameters:
3041: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3042: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3043: identity).
3045: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3046: .vb
3047: const PetscInt **perms;
3048: const PetscScalar **rots;
3049: PetscInt lOffset;
3051: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3052: for (i = 0, lOffset = 0; i < numPoints; i++) {
3053: PetscInt point = points[2*i], dof, sOffset;
3054: const PetscInt *perm = perms ? perms[i] : NULL;
3055: const PetscScalar *rot = rots ? rots[i] : NULL;
3057: PetscSectionGetDof(section,point,&dof);
3058: PetscSectionGetOffset(section,point,&sOffset);
3060: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3061: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3062: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3063: lOffset += dof;
3064: }
3065: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3066: .ve
3068: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3069: .vb
3070: const PetscInt **perms;
3071: const PetscScalar **rots;
3072: PetscInt lOffset;
3074: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3075: for (i = 0, lOffset = 0; i < numPoints; i++) {
3076: PetscInt point = points[2*i], dof, sOffset;
3077: const PetscInt *perm = perms ? perms[i] : NULL;
3078: const PetscScalar *rot = rots ? rots[i] : NULL;
3080: PetscSectionGetDof(section,point,&dof);
3081: PetscSectionGetOffset(section,point,&sOff);
3083: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3084: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3085: offset += dof;
3086: }
3087: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3088: .ve
3090: Level: developer
3092: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3093: @*/
3094: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3095: {
3096: PetscSectionSym sym;
3100: if (perms) *perms = NULL;
3101: if (rots) *rots = NULL;
3102: sym = section->sym;
3103: if (sym && (perms || rots)) {
3104: SymWorkLink link;
3106: if (sym->workin) {
3107: link = sym->workin;
3108: sym->workin = sym->workin->next;
3109: } else {
3110: PetscNewLog(sym,&link);
3111: }
3112: if (numPoints > link->numPoints) {
3113: PetscInt **perms = (PetscInt **)link->perms;
3114: PetscScalar **rots = (PetscScalar **) link->rots;
3115: PetscFree2(perms,rots);
3116: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3117: link->numPoints = numPoints;
3118: }
3119: link->next = sym->workout;
3120: sym->workout = link;
3121: PetscArrayzero((PetscInt**)link->perms,numPoints);
3122: PetscArrayzero((PetscInt**)link->rots,numPoints);
3123: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3124: if (perms) *perms = link->perms;
3125: if (rots) *rots = link->rots;
3126: }
3127: return 0;
3128: }
3130: /*@C
3131: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
3133: Not Collective
3135: Input Parameters:
3136: + section - the section
3137: . numPoints - the number of points
3138: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3139: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3140: context, see DMPlexGetConeOrientation()).
3142: Output Parameters:
3143: + perms - The permutations for the given orientations: set to NULL at conclusion
3144: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3146: Level: developer
3148: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3149: @*/
3150: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3151: {
3152: PetscSectionSym sym;
3155: sym = section->sym;
3156: if (sym && (perms || rots)) {
3157: SymWorkLink *p,link;
3159: for (p=&sym->workout; (link=*p); p=&link->next) {
3160: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3161: *p = link->next;
3162: link->next = sym->workin;
3163: sym->workin = link;
3164: if (perms) *perms = NULL;
3165: if (rots) *rots = NULL;
3166: return 0;
3167: }
3168: }
3169: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3170: }
3171: return 0;
3172: }
3174: /*@C
3175: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
3177: Not Collective
3179: Input Parameters:
3180: + section - the section
3181: . field - the field of the section
3182: . numPoints - the number of points
3183: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3184: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3185: context, see DMPlexGetConeOrientation()).
3187: Output Parameters:
3188: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3189: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3190: identity).
3192: Level: developer
3194: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3195: @*/
3196: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3197: {
3200: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3201: return 0;
3202: }
3204: /*@C
3205: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
3207: Not Collective
3209: Input Parameters:
3210: + section - the section
3211: . field - the field number
3212: . numPoints - the number of points
3213: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3214: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3215: context, see DMPlexGetConeOrientation()).
3217: Output Parameters:
3218: + perms - The permutations for the given orientations: set to NULL at conclusion
3219: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3221: Level: developer
3223: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3224: @*/
3225: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3226: {
3229: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3230: return 0;
3231: }
3233: /*@
3234: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3236: Not Collective
3238: Input Parameter:
3239: . sym - the PetscSectionSym
3241: Output Parameter:
3242: . nsym - the equivalent symmetries
3244: Level: developer
3246: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
3247: @*/
3248: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3249: {
3252: if (sym->ops->copy) (*sym->ops->copy)(sym, nsym);
3253: return 0;
3254: }
3256: /*@
3257: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input SF
3259: Collective
3261: Input Parameters:
3262: + sym - the PetscSectionSym
3263: - migrationSF - the distribution map from roots to leaves
3265: Output Parameters:
3266: . dsym - the redistributed symmetries
3268: Level: developer
3270: .seealso: PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym(), PetscSectionSymLabelSetStratum(), PetscSectionGetPointSyms()
3271: @*/
3272: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3273: {
3277: if (sym->ops->distribute) (*sym->ops->distribute)(sym, migrationSF, dsym);
3278: return 0;
3279: }
3281: /*@
3282: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3284: Not Collective
3286: Input Parameter:
3287: . s - the global PetscSection
3289: Output Parameters:
3290: . flg - the flag
3292: Level: developer
3294: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3295: @*/
3296: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3297: {
3299: *flg = s->useFieldOff;
3300: return 0;
3301: }
3303: /*@
3304: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3306: Not Collective
3308: Input Parameters:
3309: + s - the global PetscSection
3310: - flg - the flag
3312: Level: developer
3314: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3315: @*/
3316: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3317: {
3319: s->useFieldOff = flg;
3320: return 0;
3321: }
3323: #define PetscSectionExpandPoints_Loop(TYPE) \
3324: { \
3325: PetscInt i, n, o0, o1, size; \
3326: TYPE *a0 = (TYPE*)origArray, *a1; \
3327: PetscSectionGetStorageSize(s, &size); \
3328: PetscMalloc1(size, &a1); \
3329: for (i=0; i<npoints; i++) { \
3330: PetscSectionGetOffset(origSection, points_[i], &o0); \
3331: PetscSectionGetOffset(s, i, &o1); \
3332: PetscSectionGetDof(s, i, &n); \
3333: PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3334: } \
3335: *newArray = (void*)a1; \
3336: }
3338: /*@
3339: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3341: Not Collective
3343: Input Parameters:
3344: + origSection - the PetscSection describing the layout of the array
3345: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3346: . origArray - the array; its size must be equal to the storage size of origSection
3347: - points - IS with points to extract; its indices must lie in the chart of origSection
3349: Output Parameters:
3350: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3351: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3353: Level: developer
3355: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3356: @*/
3357: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3358: {
3359: PetscSection s;
3360: const PetscInt *points_;
3361: PetscInt i, n, npoints, pStart, pEnd;
3362: PetscMPIInt unitsize;
3369: MPI_Type_size(dataType, &unitsize);
3370: ISGetLocalSize(points, &npoints);
3371: ISGetIndices(points, &points_);
3372: PetscSectionGetChart(origSection, &pStart, &pEnd);
3373: PetscSectionCreate(PETSC_COMM_SELF, &s);
3374: PetscSectionSetChart(s, 0, npoints);
3375: for (i=0; i<npoints; i++) {
3377: PetscSectionGetDof(origSection, points_[i], &n);
3378: PetscSectionSetDof(s, i, n);
3379: }
3380: PetscSectionSetUp(s);
3381: if (newArray) {
3382: if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);}
3383: else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);}
3384: else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);}
3385: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3386: }
3387: if (newSection) {
3388: *newSection = s;
3389: } else {
3390: PetscSectionDestroy(&s);
3391: }
3392: ISRestoreIndices(points, &points_);
3393: return 0;
3394: }