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` 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: .vb
24: PetscSectionCreate(MPI_Comm,PetscSection *);!
25: PetscSectionSetNumFields(PetscSection, numFields);
26: PetscSectionSetChart(PetscSection,low,high);
27: PetscSectionSetDof(PetscSection,point,numdof);
28: PetscSectionSetUp(PetscSection);
29: PetscSectionGetOffset(PetscSection,point,PetscInt *);
30: PetscSectionDestroy(PetscSection);
31: .ve
33: The `PetscSection` object and methods are intended to be used in the PETSc `Vec` and `Mat` implementations.
35: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionDestroy()`
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
40: ISInitializePackage();
42: PetscHeaderCreate(*s, PETSC_SECTION_CLASSID, "PetscSection", "Section", "IS", comm, PetscSectionDestroy, PetscSectionView);
44: (*s)->pStart = -1;
45: (*s)->pEnd = -1;
46: (*s)->perm = NULL;
47: (*s)->pointMajor = PETSC_TRUE;
48: (*s)->includesConstraints = PETSC_TRUE;
49: (*s)->atlasDof = NULL;
50: (*s)->atlasOff = NULL;
51: (*s)->bc = NULL;
52: (*s)->bcIndices = NULL;
53: (*s)->setup = PETSC_FALSE;
54: (*s)->numFields = 0;
55: (*s)->fieldNames = NULL;
56: (*s)->field = NULL;
57: (*s)->useFieldOff = PETSC_FALSE;
58: (*s)->compNames = NULL;
59: (*s)->clObj = NULL;
60: (*s)->clHash = NULL;
61: (*s)->clSection = NULL;
62: (*s)->clPoints = NULL;
63: PetscSectionInvalidateMaxDof_Internal(*s);
64: return 0;
65: }
67: /*@
68: PetscSectionCopy - Creates a shallow (if possible) copy of the `PetscSection`
70: Collective
72: Input Parameter:
73: . section - the `PetscSection`
75: Output Parameter:
76: . newSection - the copy
78: Level: intermediate
80: Developer Note:
81: What exactly does shallow mean in this context?
83: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
84: @*/
85: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
86: {
87: PetscSectionSym sym;
88: IS perm;
89: PetscInt numFields, f, c, pStart, pEnd, p;
93: PetscSectionReset(newSection);
94: PetscSectionGetNumFields(section, &numFields);
95: if (numFields) PetscSectionSetNumFields(newSection, numFields);
96: for (f = 0; f < numFields; ++f) {
97: const char *fieldName = NULL, *compName = NULL;
98: PetscInt numComp = 0;
100: PetscSectionGetFieldName(section, f, &fieldName);
101: PetscSectionSetFieldName(newSection, f, fieldName);
102: PetscSectionGetFieldComponents(section, f, &numComp);
103: PetscSectionSetFieldComponents(newSection, f, numComp);
104: for (c = 0; c < numComp; ++c) {
105: PetscSectionGetComponentName(section, f, c, &compName);
106: PetscSectionSetComponentName(newSection, f, c, compName);
107: }
108: PetscSectionGetFieldSym(section, f, &sym);
109: PetscSectionSetFieldSym(newSection, f, sym);
110: }
111: PetscSectionGetChart(section, &pStart, &pEnd);
112: PetscSectionSetChart(newSection, pStart, pEnd);
113: PetscSectionGetPermutation(section, &perm);
114: PetscSectionSetPermutation(newSection, perm);
115: PetscSectionGetSym(section, &sym);
116: PetscSectionSetSym(newSection, sym);
117: for (p = pStart; p < pEnd; ++p) {
118: PetscInt dof, cdof, fcdof = 0;
120: PetscSectionGetDof(section, p, &dof);
121: PetscSectionSetDof(newSection, p, dof);
122: PetscSectionGetConstraintDof(section, p, &cdof);
123: if (cdof) PetscSectionSetConstraintDof(newSection, p, cdof);
124: for (f = 0; f < numFields; ++f) {
125: PetscSectionGetFieldDof(section, p, f, &dof);
126: PetscSectionSetFieldDof(newSection, p, f, dof);
127: if (cdof) {
128: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
129: if (fcdof) PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);
130: }
131: }
132: }
133: PetscSectionSetUp(newSection);
134: for (p = pStart; p < pEnd; ++p) {
135: PetscInt off, cdof, fcdof = 0;
136: const PetscInt *cInd;
138: /* Must set offsets in case they do not agree with the prefix sums */
139: PetscSectionGetOffset(section, p, &off);
140: PetscSectionSetOffset(newSection, p, off);
141: PetscSectionGetConstraintDof(section, p, &cdof);
142: if (cdof) {
143: PetscSectionGetConstraintIndices(section, p, &cInd);
144: PetscSectionSetConstraintIndices(newSection, p, cInd);
145: for (f = 0; f < numFields; ++f) {
146: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
147: if (fcdof) {
148: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
149: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
150: }
151: }
152: }
153: }
154: return 0;
155: }
157: /*@
158: PetscSectionClone - Creates a shallow (if possible) copy of the `PetscSection`
160: Collective
162: Input Parameter:
163: . section - the `PetscSection`
165: Output Parameter:
166: . newSection - the copy
168: Level: beginner
170: Developer Note:
171: With standard PETSc terminology this should be called `PetscSectionDuplicate()`
173: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
174: @*/
175: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
176: {
179: PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection);
180: PetscSectionCopy(section, *newSection);
181: return 0;
182: }
184: /*@
185: PetscSectionSetFromOptions - sets parameters in a `PetscSection` from the options database
187: Collective
189: Input Parameter:
190: . section - the `PetscSection`
192: Options Database:
193: . -petscsection_point_major - `PETSC_TRUE` for point-major order
195: Level: intermediate
197: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
198: @*/
199: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
200: {
202: PetscObjectOptionsBegin((PetscObject)s);
203: PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
204: /* process any options handlers added with PetscObjectAddOptionsHandler() */
205: PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject);
206: PetscOptionsEnd();
207: PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view");
208: return 0;
209: }
211: /*@
212: PetscSectionCompare - Compares two sections
214: Collective
216: Input Parameters:
217: + s1 - the first `PetscSection`
218: - s2 - the second `PetscSection`
220: Output Parameter:
221: . congruent - `PETSC_TRUE` if the two sections are congruent, `PETSC_FALSE` otherwise
223: Level: intermediate
225: Note:
226: Field names are disregarded.
228: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
229: @*/
230: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
231: {
232: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
233: const PetscInt *idx1, *idx2;
234: IS perm1, perm2;
235: PetscBool flg;
236: PetscMPIInt mflg;
241: flg = PETSC_FALSE;
243: MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg);
244: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
245: *congruent = PETSC_FALSE;
246: return 0;
247: }
249: PetscSectionGetChart(s1, &pStart, &pEnd);
250: PetscSectionGetChart(s2, &n1, &n2);
251: if (pStart != n1 || pEnd != n2) goto not_congruent;
253: PetscSectionGetPermutation(s1, &perm1);
254: PetscSectionGetPermutation(s2, &perm2);
255: if (perm1 && perm2) {
256: ISEqual(perm1, perm2, congruent);
257: if (!(*congruent)) goto not_congruent;
258: } else if (perm1 != perm2) goto not_congruent;
260: for (p = pStart; p < pEnd; ++p) {
261: PetscSectionGetOffset(s1, p, &n1);
262: PetscSectionGetOffset(s2, p, &n2);
263: if (n1 != n2) goto not_congruent;
265: PetscSectionGetDof(s1, p, &n1);
266: PetscSectionGetDof(s2, p, &n2);
267: if (n1 != n2) goto not_congruent;
269: PetscSectionGetConstraintDof(s1, p, &ncdof);
270: PetscSectionGetConstraintDof(s2, p, &n2);
271: if (ncdof != n2) goto not_congruent;
273: PetscSectionGetConstraintIndices(s1, p, &idx1);
274: PetscSectionGetConstraintIndices(s2, p, &idx2);
275: PetscArraycmp(idx1, idx2, ncdof, congruent);
276: if (!(*congruent)) goto not_congruent;
277: }
279: PetscSectionGetNumFields(s1, &nfields);
280: PetscSectionGetNumFields(s2, &n2);
281: if (nfields != n2) goto not_congruent;
283: for (f = 0; f < nfields; ++f) {
284: PetscSectionGetFieldComponents(s1, f, &n1);
285: PetscSectionGetFieldComponents(s2, f, &n2);
286: if (n1 != n2) goto not_congruent;
288: for (p = pStart; p < pEnd; ++p) {
289: PetscSectionGetFieldOffset(s1, p, f, &n1);
290: PetscSectionGetFieldOffset(s2, p, f, &n2);
291: if (n1 != n2) goto not_congruent;
293: PetscSectionGetFieldDof(s1, p, f, &n1);
294: PetscSectionGetFieldDof(s2, p, f, &n2);
295: if (n1 != n2) goto not_congruent;
297: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
298: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
299: if (nfcdof != n2) goto not_congruent;
301: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
302: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
303: PetscArraycmp(idx1, idx2, nfcdof, congruent);
304: if (!(*congruent)) goto not_congruent;
305: }
306: }
308: flg = PETSC_TRUE;
309: not_congruent:
310: MPIU_Allreduce(&flg, congruent, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1));
311: return 0;
312: }
314: /*@
315: PetscSectionGetNumFields - Returns the number of fields in a `PetscSection`, or 0 if no fields were defined.
317: Not Collective
319: Input Parameter:
320: . s - the `PetscSection`
322: Output Parameter:
323: . numFields - the number of fields defined, or 0 if none were defined
325: Level: intermediate
327: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
328: @*/
329: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
330: {
333: *numFields = s->numFields;
334: return 0;
335: }
337: /*@
338: PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`
340: Not Collective
342: Input Parameters:
343: + s - the PetscSection
344: - numFields - the number of fields
346: Level: intermediate
348: .seealso: [PetscSection](sec_petscsection), `PetscSection()`, `PetscSectionGetNumFields()`
349: @*/
350: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
351: {
352: PetscInt f;
356: PetscSectionReset(s);
358: s->numFields = numFields;
359: PetscMalloc1(s->numFields, &s->numFieldComponents);
360: PetscMalloc1(s->numFields, &s->fieldNames);
361: PetscMalloc1(s->numFields, &s->compNames);
362: PetscMalloc1(s->numFields, &s->field);
363: for (f = 0; f < s->numFields; ++f) {
364: char name[64];
366: s->numFieldComponents[f] = 1;
368: PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]);
369: PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f);
370: PetscStrallocpy(name, (char **)&s->fieldNames[f]);
371: PetscSNPrintf(name, 64, "Component_0");
372: PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
373: PetscStrallocpy(name, (char **)&s->compNames[f][0]);
374: }
375: return 0;
376: }
378: /*@C
379: PetscSectionGetFieldName - Returns the name of a field in the `PetscSection`
381: Not Collective
383: Input Parameters:
384: + s - the `PetscSection`
385: - field - the field number
387: Output Parameter:
388: . fieldName - the field name
390: Level: intermediate
392: Note:
393: Will error if the field number is out of range
395: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
396: @*/
397: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
398: {
401: PetscSectionCheckValidField(field, s->numFields);
402: *fieldName = s->fieldNames[field];
403: return 0;
404: }
406: /*@C
407: PetscSectionSetFieldName - Sets the name of a field in the `PetscSection`
409: Not Collective
411: Input Parameters:
412: + s - the `PetscSection`
413: . field - the field number
414: - fieldName - the field name
416: Level: intermediate
418: Note:
419: Will error if the field number is out of range
421: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
422: @*/
423: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
424: {
427: PetscSectionCheckValidField(field, s->numFields);
428: PetscFree(s->fieldNames[field]);
429: PetscStrallocpy(fieldName, (char **)&s->fieldNames[field]);
430: return 0;
431: }
433: /*@C
434: PetscSectionGetComponentName - Gets the name of a field component in the `PetscSection`
436: Not Collective
438: Input Parameters:
439: + s - the `PetscSection`
440: . field - the field number
441: - comp - the component number
443: Ouput Parameter:
444: . compName - the component name
446: Level: intermediate
448: Note:
449: Will error if the field or component number do not exist
451: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
452: `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
453: @*/
454: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
455: {
458: PetscSectionCheckValidField(field, s->numFields);
459: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
460: *compName = s->compNames[field][comp];
461: return 0;
462: }
464: /*@C
465: PetscSectionSetComponentName - Sets the name of a field component in the `PetscSection`
467: Not Collective
469: Input Parameters:
470: + s - the PetscSection
471: . field - the field number
472: . comp - the component number
473: - compName - the component name
475: Level: intermediate
477: Note:
478: Will error if the field or component number do not exist
480: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
481: `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
482: @*/
483: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
484: {
487: PetscSectionCheckValidField(field, s->numFields);
488: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
489: PetscFree(s->compNames[field][comp]);
490: PetscStrallocpy(compName, (char **)&s->compNames[field][comp]);
491: return 0;
492: }
494: /*@
495: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
497: Not Collective
499: Input Parameters:
500: + s - the `PetscSection`
501: - field - the field number
503: Output Parameter:
504: . numComp - the number of field components
506: Level: intermediate
508: Developer Note:
509: This function is misnamed. There is a Num in `PetscSectionGetNumFields()` but not in this name
511: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`
512: @*/
513: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
514: {
517: PetscSectionCheckValidField(field, s->numFields);
518: *numComp = s->numFieldComponents[field];
519: return 0;
520: }
522: /*@
523: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
525: Not Collective
527: Input Parameters:
528: + s - the PetscSection
529: . field - the field number
530: - numComp - the number of field components
532: Level: intermediate
534: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionGetNumFields()`
535: @*/
536: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
537: {
538: PetscInt c;
541: PetscSectionCheckValidField(field, s->numFields);
542: if (s->compNames) {
543: for (c = 0; c < s->numFieldComponents[field]; ++c) PetscFree(s->compNames[field][c]);
544: PetscFree(s->compNames[field]);
545: }
547: s->numFieldComponents[field] = numComp;
548: if (numComp) {
549: PetscMalloc1(numComp, (char ***)&s->compNames[field]);
550: for (c = 0; c < numComp; ++c) {
551: char name[64];
553: PetscSNPrintf(name, 64, "%" PetscInt_FMT, c);
554: PetscStrallocpy(name, (char **)&s->compNames[field][c]);
555: }
556: }
557: return 0;
558: }
560: /*@
561: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points (indices) lie for this `PetscSection`
563: Not Collective
565: Input Parameter:
566: . s - the `PetscSection`
568: Output Parameters:
569: + pStart - the first point
570: - pEnd - one past the last point
572: Level: intermediate
574: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
575: @*/
576: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
577: {
579: if (pStart) *pStart = s->pStart;
580: if (pEnd) *pEnd = s->pEnd;
581: return 0;
582: }
584: /*@
585: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points (indices) lie for this `PetscSection`
587: Not Collective
589: Input Parameters:
590: + s - the PetscSection
591: . pStart - the first point
592: - pEnd - one past the last point
594: Level: intermediate
596: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`
597: @*/
598: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
599: {
600: PetscInt f;
603: if (pStart == s->pStart && pEnd == s->pEnd) return 0;
604: /* Cannot Reset() because it destroys field information */
605: s->setup = PETSC_FALSE;
606: PetscSectionDestroy(&s->bc);
607: PetscFree(s->bcIndices);
608: PetscFree2(s->atlasDof, s->atlasOff);
610: s->pStart = pStart;
611: s->pEnd = pEnd;
612: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
613: PetscArrayzero(s->atlasDof, pEnd - pStart);
614: for (f = 0; f < s->numFields; ++f) PetscSectionSetChart(s->field[f], pStart, pEnd);
615: return 0;
616: }
618: /*@
619: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL that was set with `PetscSectionSetPermutation()`
621: Not Collective
623: Input Parameter:
624: . s - the `PetscSection`
626: Output Parameters:
627: . perm - The permutation as an `IS`
629: Level: intermediate
631: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
632: @*/
633: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
634: {
636: if (perm) {
638: *perm = s->perm;
639: }
640: return 0;
641: }
643: /*@
644: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
646: Not Collective
648: Input Parameters:
649: + s - the PetscSection
650: - perm - the permutation of points
652: Level: intermediate
654: Developer Note:
655: What purpose does this permutation serve?
657: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionGetPermutation()`, `PetscSectionCreate()`
658: @*/
659: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
660: {
664: if (s->perm != perm) {
665: ISDestroy(&s->perm);
666: if (perm) {
667: s->perm = perm;
668: PetscObjectReference((PetscObject)s->perm);
669: }
670: }
671: return 0;
672: }
674: /*@
675: PetscSectionGetPointMajor - Returns the flag for dof ordering, `PETSC_TRUE` if it is point major, `PETSC_FALSE` if it is field major
677: Not Collective
679: Input Parameter:
680: . s - the `PetscSection`
682: Output Parameter:
683: . pm - the flag for point major ordering
685: Level: intermediate
687: .seealso: [PetscSection](sec_petscsection), `PetscSection`, PetscSectionSetPointMajor()`
688: @*/
689: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
690: {
693: *pm = s->pointMajor;
694: return 0;
695: }
697: /*@
698: PetscSectionSetPointMajor - Sets the flag for dof ordering, `PETSC_TRUE` for point major, otherwise it will be field major
700: Not Collective
702: Input Parameters:
703: + s - the `PetscSection`
704: - pm - the flag for point major ordering
706: Level: intermediate
708: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`
709: @*/
710: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
711: {
714: s->pointMajor = pm;
715: return 0;
716: }
718: /*@
719: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets in the `PetscSection`.
720: The value is set with `PetscSectionSetIncludesConstraints()`
722: Not Collective
724: Input Parameter:
725: . s - the `PetscSection`
727: Output Parameter:
728: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
730: Level: intermediate
732: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
733: @*/
734: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
735: {
738: *includesConstraints = s->includesConstraints;
739: return 0;
740: }
742: /*@
743: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
745: Not Collective
747: Input Parameters:
748: + s - the PetscSection
749: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
751: Level: intermediate
753: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
754: @*/
755: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
756: {
759: s->includesConstraints = includesConstraints;
760: return 0;
761: }
763: /*@
764: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
766: Not Collective
768: Input Parameters:
769: + s - the `PetscSection`
770: - point - the point
772: Output Parameter:
773: . numDof - the number of dof
775: Level: intermediate
777: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
778: @*/
779: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
780: {
784: 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);
785: *numDof = s->atlasDof[point - s->pStart];
786: return 0;
787: }
789: /*@
790: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
792: Not Collective
794: Input Parameters:
795: + s - the `PetscSection`
796: . point - the point
797: - numDof - the number of dof
799: Level: intermediate
801: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
802: @*/
803: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
804: {
806: 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);
807: s->atlasDof[point - s->pStart] = numDof;
808: PetscSectionInvalidateMaxDof_Internal(s);
809: return 0;
810: }
812: /*@
813: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
815: Not Collective
817: Input Parameters:
818: + s - the `PetscSection`
819: . point - the point
820: - numDof - the number of additional dof
822: Level: intermediate
824: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
825: @*/
826: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
827: {
830: 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);
831: s->atlasDof[point - s->pStart] += numDof;
832: PetscSectionInvalidateMaxDof_Internal(s);
833: return 0;
834: }
836: /*@
837: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
839: Not Collective
841: Input Parameters:
842: + s - the `PetscSection`
843: . point - the point
844: - field - the field
846: Output Parameter:
847: . numDof - the number of dof
849: Level: intermediate
851: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
852: @*/
853: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
854: {
857: PetscSectionCheckValidField(field, s->numFields);
858: PetscSectionGetDof(s->field[field], point, numDof);
859: return 0;
860: }
862: /*@
863: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
865: Not Collective
867: Input Parameters:
868: + s - the `PetscSection`
869: . point - the point
870: . field - the field
871: - numDof - the number of dof
873: Level: intermediate
875: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
876: @*/
877: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
878: {
880: PetscSectionCheckValidField(field, s->numFields);
881: PetscSectionSetDof(s->field[field], point, numDof);
882: return 0;
883: }
885: /*@
886: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
888: Not Collective
890: Input Parameters:
891: + s - the `PetscSection`
892: . point - the point
893: . field - the field
894: - numDof - the number of dof
896: Level: intermediate
898: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
899: @*/
900: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
901: {
903: PetscSectionCheckValidField(field, s->numFields);
904: PetscSectionAddDof(s->field[field], point, numDof);
905: return 0;
906: }
908: /*@
909: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
911: Not Collective
913: Input Parameters:
914: + s - the `PetscSection`
915: - point - the point
917: Output Parameter:
918: . numDof - the number of dof which are fixed by constraints
920: Level: intermediate
922: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
923: @*/
924: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
925: {
928: if (s->bc) {
929: PetscSectionGetDof(s->bc, point, numDof);
930: } else *numDof = 0;
931: return 0;
932: }
934: /*@
935: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
937: Not Collective
939: Input Parameters:
940: + s - the `PetscSection`
941: . point - the point
942: - numDof - the number of dof which are fixed by constraints
944: Level: intermediate
946: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
947: @*/
948: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
949: {
951: if (numDof) {
952: PetscSectionCheckConstraints_Private(s);
953: PetscSectionSetDof(s->bc, point, numDof);
954: }
955: return 0;
956: }
958: /*@
959: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
961: Not Collective
963: Input Parameters:
964: + s - the `PetscSection`
965: . point - the point
966: - numDof - the number of additional dof which are fixed by constraints
968: Level: intermediate
970: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
971: @*/
972: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
973: {
975: if (numDof) {
976: PetscSectionCheckConstraints_Private(s);
977: PetscSectionAddDof(s->bc, point, numDof);
978: }
979: return 0;
980: }
982: /*@
983: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
985: Not Collective
987: Input Parameters:
988: + s - the `PetscSection`
989: . point - the point
990: - field - the field
992: Output Parameter:
993: . numDof - the number of dof which are fixed by constraints
995: Level: intermediate
997: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
998: @*/
999: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1000: {
1003: PetscSectionCheckValidField(field, s->numFields);
1004: PetscSectionGetConstraintDof(s->field[field], point, numDof);
1005: return 0;
1006: }
1008: /*@
1009: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1011: Not Collective
1013: Input Parameters:
1014: + s - the `PetscSection`
1015: . point - the point
1016: . field - the field
1017: - numDof - the number of dof which are fixed by constraints
1019: Level: intermediate
1021: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1022: @*/
1023: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1024: {
1026: PetscSectionCheckValidField(field, s->numFields);
1027: PetscSectionSetConstraintDof(s->field[field], point, numDof);
1028: return 0;
1029: }
1031: /*@
1032: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1034: Not Collective
1036: Input Parameters:
1037: + s - the `PetscSection`
1038: . point - the point
1039: . field - the field
1040: - numDof - the number of additional dof which are fixed by constraints
1042: Level: intermediate
1044: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1045: @*/
1046: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1047: {
1049: PetscSectionCheckValidField(field, s->numFields);
1050: PetscSectionAddConstraintDof(s->field[field], point, numDof);
1051: return 0;
1052: }
1054: /*@
1055: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1057: Not Collective
1059: Input Parameter:
1060: . s - the `PetscSection`
1062: Level: advanced
1064: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1065: @*/
1066: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1067: {
1069: if (s->bc) {
1070: const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;
1072: PetscSectionSetUp(s->bc);
1073: PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices);
1074: }
1075: return 0;
1076: }
1078: /*@
1079: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
1081: Not Collective
1083: Input Parameter:
1084: . s - the `PetscSection`
1086: Level: intermediate
1088: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1089: @*/
1090: PetscErrorCode PetscSectionSetUp(PetscSection s)
1091: {
1092: const PetscInt *pind = NULL;
1093: PetscInt offset = 0, foff, p, f;
1096: if (s->setup) return 0;
1097: s->setup = PETSC_TRUE;
1098: /* Set offsets and field offsets for all points */
1099: /* Assume that all fields have the same chart */
1101: if (s->perm) ISGetIndices(s->perm, &pind);
1102: if (s->pointMajor) {
1103: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1104: const PetscInt q = pind ? pind[p] : p;
1106: /* Set point offset */
1107: s->atlasOff[q] = offset;
1108: offset += s->atlasDof[q];
1109: /* Set field offset */
1110: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1111: PetscSection sf = s->field[f];
1113: sf->atlasOff[q] = foff;
1114: foff += sf->atlasDof[q];
1115: }
1116: }
1117: } else {
1118: /* Set field offsets for all points */
1119: for (f = 0; f < s->numFields; ++f) {
1120: PetscSection sf = s->field[f];
1122: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1123: const PetscInt q = pind ? pind[p] : p;
1125: sf->atlasOff[q] = offset;
1126: offset += sf->atlasDof[q];
1127: }
1128: }
1129: /* Disable point offsets since these are unused */
1130: for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1131: }
1132: if (s->perm) ISRestoreIndices(s->perm, &pind);
1133: /* Setup BC sections */
1134: PetscSectionSetUpBC(s);
1135: for (f = 0; f < s->numFields; ++f) PetscSectionSetUpBC(s->field[f]);
1136: return 0;
1137: }
1139: /*@
1140: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1142: Not Collective
1144: Input Parameters:
1145: . s - the `PetscSection`
1147: Output Parameter:
1148: . maxDof - the maximum dof
1150: Level: intermediate
1152: Note:
1153: The returned number is up-to-date without need for `PetscSectionSetUp()`.
1155: Developer Note:
1156: The returned number is calculated lazily and stashed.
1158: A call to `PetscSectionInvalidateMaxDof_Internal()` invalidates the stashed value.
1160: `PetscSectionInvalidateMaxDof_Internal()` is called in `PetscSectionSetDof()`, `PetscSectionAddDof()` and `PetscSectionReset()`
1162: It should also be called every time atlasDof is modified directly.
1164: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1165: @*/
1166: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1167: {
1168: PetscInt p;
1172: if (s->maxDof == PETSC_MIN_INT) {
1173: s->maxDof = 0;
1174: for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1175: }
1176: *maxDof = s->maxDof;
1177: return 0;
1178: }
1180: /*@
1181: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the 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 the dofs
1191: Level: intermediate
1193: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1194: @*/
1195: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1196: {
1197: PetscInt p, n = 0;
1201: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1202: *size = n;
1203: return 0;
1204: }
1206: /*@
1207: PetscSectionGetConstrainedStorageSize - Return the size of an array or local `Vec` capable of holding all unconstrained degrees of freedom.
1209: Not Collective
1211: Input Parameter:
1212: . s - the `PetscSection`
1214: Output Parameter:
1215: . size - the size of an array which can hold all unconstrained dofs
1217: Level: intermediate
1219: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1220: @*/
1221: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1222: {
1223: PetscInt p, n = 0;
1227: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1228: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1229: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1230: }
1231: *size = n;
1232: return 0;
1233: }
1235: /*@
1236: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1237: the local section and a `PetscSF` describing the section point overlap.
1239: Input Parameters:
1240: + s - The `PetscSection` for the local field layout
1241: . sf - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1242: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1243: - localOffsets - If `PETSC_TRUE`, use local rather than global offsets for the points
1245: Output Parameter:
1246: . gsection - The `PetscSection` for the global field layout
1248: Level: intermediate
1250: Note:
1251: This gives negative sizes and offsets to points not owned by this process
1253: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1254: @*/
1255: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1256: {
1257: PetscSection gs;
1258: const PetscInt *pind = NULL;
1259: PetscInt *recv = NULL, *neg = NULL;
1260: PetscInt pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1261: PetscInt numFields, f, numComponents;
1269: PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs);
1270: PetscSectionGetNumFields(s, &numFields);
1271: if (numFields > 0) PetscSectionSetNumFields(gs, numFields);
1272: PetscSectionGetChart(s, &pStart, &pEnd);
1273: PetscSectionSetChart(gs, pStart, pEnd);
1274: gs->includesConstraints = includeConstraints;
1275: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1276: nlocal = nroots; /* The local/leaf space matches global/root space */
1277: /* Must allocate for all points visible to SF, which may be more than this section */
1278: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1279: PetscSFGetLeafRange(sf, NULL, &maxleaf);
1282: PetscMalloc2(nroots, &neg, nlocal, &recv);
1283: PetscArrayzero(neg, nroots);
1284: }
1285: /* Mark all local points with negative dof */
1286: for (p = pStart; p < pEnd; ++p) {
1287: PetscSectionGetDof(s, p, &dof);
1288: PetscSectionSetDof(gs, p, dof);
1289: PetscSectionGetConstraintDof(s, p, &cdof);
1290: if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(gs, p, cdof);
1291: if (neg) neg[p] = -(dof + 1);
1292: }
1293: PetscSectionSetUpBC(gs);
1294: 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]);
1295: if (nroots >= 0) {
1296: PetscArrayzero(recv, nlocal);
1297: PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1298: PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1299: for (p = pStart; p < pEnd; ++p) {
1300: if (recv[p] < 0) {
1301: gs->atlasDof[p - pStart] = recv[p];
1302: PetscSectionGetDof(s, p, &dof);
1304: }
1305: }
1306: }
1307: /* Calculate new sizes, get process offset, and calculate point offsets */
1308: if (s->perm) ISGetIndices(s->perm, &pind);
1309: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1310: const PetscInt q = pind ? pind[p] : p;
1312: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1313: gs->atlasOff[q] = off;
1314: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1315: }
1316: if (!localOffsets) {
1317: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
1318: globalOff -= off;
1319: }
1320: for (p = pStart, off = 0; p < pEnd; ++p) {
1321: gs->atlasOff[p - pStart] += globalOff;
1322: if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1323: }
1324: if (s->perm) ISRestoreIndices(s->perm, &pind);
1325: /* Put in negative offsets for ghost points */
1326: if (nroots >= 0) {
1327: PetscArrayzero(recv, nlocal);
1328: PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1329: PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1330: for (p = pStart; p < pEnd; ++p) {
1331: if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1332: }
1333: }
1334: PetscFree2(neg, recv);
1335: /* Set field dofs/offsets/constraints */
1336: for (f = 0; f < numFields; ++f) {
1337: gs->field[f]->includesConstraints = includeConstraints;
1338: PetscSectionGetFieldComponents(s, f, &numComponents);
1339: PetscSectionSetFieldComponents(gs, f, numComponents);
1340: }
1341: for (p = pStart; p < pEnd; ++p) {
1342: PetscSectionGetOffset(gs, p, &off);
1343: for (f = 0, foff = off; f < numFields; ++f) {
1344: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1345: if (!includeConstraints && cdof > 0) PetscSectionSetFieldConstraintDof(gs, p, f, cdof);
1346: PetscSectionGetFieldDof(s, p, f, &dof);
1347: PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof);
1348: PetscSectionSetFieldOffset(gs, p, f, foff);
1349: PetscSectionGetFieldConstraintDof(gs, p, f, &cdof);
1350: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1351: }
1352: }
1353: for (f = 0; f < numFields; ++f) {
1354: PetscSection gfs = gs->field[f];
1356: PetscSectionSetUpBC(gfs);
1357: 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]);
1358: }
1359: gs->setup = PETSC_TRUE;
1360: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1361: *gsection = gs;
1362: return 0;
1363: }
1365: /*@
1366: PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the global field layout using
1367: the local section and an `PetscSF` describing the section point overlap.
1369: Input Parameters:
1370: + s - The `PetscSection` for the local field layout
1371: . sf - The `PetscSF` describing parallel layout of the section points
1372: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1373: . numExcludes - The number of exclusion ranges
1374: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1376: Output Parameter:
1377: . gsection - The `PetscSection` for the global field layout
1379: Level: advanced
1381: Note:
1382: This gives negative sizes and offsets to points not owned by this process
1384: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1385: @*/
1386: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1387: {
1388: const PetscInt *pind = NULL;
1389: PetscInt *neg = NULL, *tmpOff = NULL;
1390: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1395: PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection);
1396: PetscSectionGetChart(s, &pStart, &pEnd);
1397: PetscSectionSetChart(*gsection, pStart, pEnd);
1398: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1399: if (nroots >= 0) {
1401: PetscCalloc1(nroots, &neg);
1402: if (nroots > pEnd - pStart) {
1403: PetscCalloc1(nroots, &tmpOff);
1404: } else {
1405: tmpOff = &(*gsection)->atlasDof[-pStart];
1406: }
1407: }
1408: /* Mark ghost points with negative dof */
1409: for (p = pStart; p < pEnd; ++p) {
1410: for (e = 0; e < numExcludes; ++e) {
1411: if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1412: PetscSectionSetDof(*gsection, p, 0);
1413: break;
1414: }
1415: }
1416: if (e < numExcludes) continue;
1417: PetscSectionGetDof(s, p, &dof);
1418: PetscSectionSetDof(*gsection, p, dof);
1419: PetscSectionGetConstraintDof(s, p, &cdof);
1420: if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(*gsection, p, cdof);
1421: if (neg) neg[p] = -(dof + 1);
1422: }
1423: PetscSectionSetUpBC(*gsection);
1424: if (nroots >= 0) {
1425: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1426: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1427: if (nroots > pEnd - pStart) {
1428: for (p = pStart; p < pEnd; ++p) {
1429: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1430: }
1431: }
1432: }
1433: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1434: if (s->perm) ISGetIndices(s->perm, &pind);
1435: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1436: const PetscInt q = pind ? pind[p] : p;
1438: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1439: (*gsection)->atlasOff[q] = off;
1440: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1441: }
1442: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
1443: globalOff -= off;
1444: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1445: (*gsection)->atlasOff[p] += globalOff;
1446: if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1447: }
1448: if (s->perm) ISRestoreIndices(s->perm, &pind);
1449: /* Put in negative offsets for ghost points */
1450: if (nroots >= 0) {
1451: if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1452: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1453: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1454: if (nroots > pEnd - pStart) {
1455: for (p = pStart; p < pEnd; ++p) {
1456: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1457: }
1458: }
1459: }
1460: if (nroots >= 0 && nroots > pEnd - pStart) PetscFree(tmpOff);
1461: PetscFree(neg);
1462: return 0;
1463: }
1465: /*@
1466: PetscSectionGetPointLayout - Get the `PetscLayout` associated with the section points
1468: Collective
1470: Input Parameters:
1471: + comm - The `MPI_Comm`
1472: - s - The `PetscSection`
1474: Output Parameter:
1475: . layout - The point layout for the section
1477: Level: advanced
1479: Note:
1480: This is usually called for the default global section.
1482: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1483: @*/
1484: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1485: {
1486: PetscInt pStart, pEnd, p, localSize = 0;
1488: PetscSectionGetChart(s, &pStart, &pEnd);
1489: for (p = pStart; p < pEnd; ++p) {
1490: PetscInt dof;
1492: PetscSectionGetDof(s, p, &dof);
1493: if (dof >= 0) ++localSize;
1494: }
1495: PetscLayoutCreate(comm, layout);
1496: PetscLayoutSetLocalSize(*layout, localSize);
1497: PetscLayoutSetBlockSize(*layout, 1);
1498: PetscLayoutSetUp(*layout);
1499: return 0;
1500: }
1502: /*@
1503: PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs.
1505: Collective
1507: Input Parameters:
1508: + comm - The `MPI_Comm`
1509: - s - The `PetscSection`
1511: Output Parameter:
1512: . layout - The dof layout for the section
1514: Level: advanced
1516: Note:
1517: This is usually called for the default global section.
1519: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1520: @*/
1521: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1522: {
1523: PetscInt pStart, pEnd, p, localSize = 0;
1527: PetscSectionGetChart(s, &pStart, &pEnd);
1528: for (p = pStart; p < pEnd; ++p) {
1529: PetscInt dof, cdof;
1531: PetscSectionGetDof(s, p, &dof);
1532: PetscSectionGetConstraintDof(s, p, &cdof);
1533: if (dof - cdof > 0) localSize += dof - cdof;
1534: }
1535: PetscLayoutCreate(comm, layout);
1536: PetscLayoutSetLocalSize(*layout, localSize);
1537: PetscLayoutSetBlockSize(*layout, 1);
1538: PetscLayoutSetUp(*layout);
1539: return 0;
1540: }
1542: /*@
1543: PetscSectionGetOffset - Return the offset into an array or local `Vec` for the dof associated with the given point.
1545: Not Collective
1547: Input Parameters:
1548: + s - the `PetscSection`
1549: - point - the point
1551: Output Parameter:
1552: . offset - the offset
1554: Level: intermediate
1556: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`
1557: @*/
1558: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1559: {
1563: *offset = s->atlasOff[point - s->pStart];
1564: return 0;
1565: }
1567: /*@
1568: PetscSectionSetOffset - Set the offset into an array or local `Vec` for the dof associated with the given point.
1570: Not Collective
1572: Input Parameters:
1573: + s - the `PetscSection`
1574: . point - the point
1575: - offset - the offset
1577: Level: intermediate
1579: Note:
1580: The user usually does not call this function, but uses `PetscSectionSetUp()`
1582: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1583: @*/
1584: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1585: {
1588: s->atlasOff[point - s->pStart] = offset;
1589: return 0;
1590: }
1592: /*@
1593: PetscSectionGetFieldOffset - Return the offset into an array or local `Vec` for the dof associated with the given point.
1595: Not Collective
1597: Input Parameters:
1598: + s - the `PetscSection`
1599: . point - the point
1600: - field - the field
1602: Output Parameter:
1603: . offset - the offset
1605: Level: intermediate
1607: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1608: @*/
1609: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1610: {
1613: PetscSectionCheckValidField(field, s->numFields);
1614: PetscSectionGetOffset(s->field[field], point, offset);
1615: return 0;
1616: }
1618: /*@
1619: PetscSectionSetFieldOffset - Set the offset into an array or local `Vec` for the dof associated with the given field at a point.
1621: Not Collective
1623: Input Parameters:
1624: + s - the PetscSection
1625: . point - the point
1626: . field - the field
1627: - offset - the offset
1629: Level: developer
1631: Note:
1632: The user usually does not call this function, but uses `PetscSectionSetUp()`
1634: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1635: @*/
1636: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1637: {
1639: PetscSectionCheckValidField(field, s->numFields);
1640: PetscSectionSetOffset(s->field[field], point, offset);
1641: return 0;
1642: }
1644: /*@
1645: PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.
1647: Not Collective
1649: Input Parameters:
1650: + s - the `PetscSection`
1651: . point - the point
1652: - field - the field
1654: Output Parameter:
1655: . offset - the offset
1657: Level: advanced
1659: Note:
1660: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1661: this point, what number is the first dof with this field.
1663: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1664: @*/
1665: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1666: {
1667: PetscInt off, foff;
1671: PetscSectionCheckValidField(field, s->numFields);
1672: PetscSectionGetOffset(s, point, &off);
1673: PetscSectionGetOffset(s->field[field], point, &foff);
1674: *offset = foff - off;
1675: return 0;
1676: }
1678: /*@
1679: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1681: Not Collective
1683: Input Parameter:
1684: . s - the `PetscSection`
1686: Output Parameters:
1687: + start - the minimum offset
1688: - end - one more than the maximum offset
1690: Level: intermediate
1692: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1693: @*/
1694: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1695: {
1696: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1699: if (s->atlasOff) {
1700: os = s->atlasOff[0];
1701: oe = s->atlasOff[0];
1702: }
1703: PetscSectionGetChart(s, &pStart, &pEnd);
1704: for (p = 0; p < pEnd - pStart; ++p) {
1705: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1707: if (off >= 0) {
1708: os = PetscMin(os, off);
1709: oe = PetscMax(oe, off + dof);
1710: }
1711: }
1712: if (start) *start = os;
1713: if (end) *end = oe;
1714: return 0;
1715: }
1717: /*@
1718: PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields
1720: Collective
1722: Input Parameters:
1723: + s - the `PetscSection`
1724: . len - the number of subfields
1725: - fields - the subfield numbers
1727: Output Parameter:
1728: . subs - the subsection
1730: Level: advanced
1732: Notes:
1733: The section offsets now refer to a new, smaller vector.
1735: This will error if a fieldnumber is out of range
1737: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
1738: @*/
1739: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1740: {
1741: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
1743: if (!len) return 0;
1747: PetscSectionGetNumFields(s, &nF);
1749: PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1750: PetscSectionSetNumFields(*subs, len);
1751: for (f = 0; f < len; ++f) {
1752: const char *name = NULL;
1753: PetscInt numComp = 0;
1755: PetscSectionGetFieldName(s, fields[f], &name);
1756: PetscSectionSetFieldName(*subs, f, name);
1757: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1758: PetscSectionSetFieldComponents(*subs, f, numComp);
1759: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1760: PetscSectionGetComponentName(s, fields[f], c, &name);
1761: PetscSectionSetComponentName(*subs, f, c, name);
1762: }
1763: }
1764: PetscSectionGetChart(s, &pStart, &pEnd);
1765: PetscSectionSetChart(*subs, pStart, pEnd);
1766: for (p = pStart; p < pEnd; ++p) {
1767: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1769: for (f = 0; f < len; ++f) {
1770: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1771: PetscSectionSetFieldDof(*subs, p, f, fdof);
1772: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1773: if (cfdof) PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);
1774: dof += fdof;
1775: cdof += cfdof;
1776: }
1777: PetscSectionSetDof(*subs, p, dof);
1778: if (cdof) PetscSectionSetConstraintDof(*subs, p, cdof);
1779: maxCdof = PetscMax(cdof, maxCdof);
1780: }
1781: PetscSectionSetUp(*subs);
1782: if (maxCdof) {
1783: PetscInt *indices;
1785: PetscMalloc1(maxCdof, &indices);
1786: for (p = pStart; p < pEnd; ++p) {
1787: PetscInt cdof;
1789: PetscSectionGetConstraintDof(*subs, p, &cdof);
1790: if (cdof) {
1791: const PetscInt *oldIndices = NULL;
1792: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1794: for (f = 0; f < len; ++f) {
1795: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1796: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1797: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1798: PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices);
1799: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
1800: numConst += cfdof;
1801: fOff += fdof;
1802: }
1803: PetscSectionSetConstraintIndices(*subs, p, indices);
1804: }
1805: }
1806: PetscFree(indices);
1807: }
1808: return 0;
1809: }
1811: /*@
1812: PetscSectionCreateSupersection - Create a new, larger section composed of the input sections
1814: Collective
1816: Input Parameters:
1817: + s - the input sections
1818: - len - the number of input sections
1820: Output Parameter:
1821: . supers - the supersection
1823: Level: advanced
1825: Note:
1826: The section offsets now refer to a new, larger vector.
1828: Developer Note:
1829: Needs to explain how the sections are composed
1831: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
1832: @*/
1833: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1834: {
1835: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1837: if (!len) return 0;
1838: for (i = 0; i < len; ++i) {
1839: PetscInt nf, pStarti, pEndi;
1841: PetscSectionGetNumFields(s[i], &nf);
1842: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1843: pStart = PetscMin(pStart, pStarti);
1844: pEnd = PetscMax(pEnd, pEndi);
1845: Nf += nf;
1846: }
1847: PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers);
1848: PetscSectionSetNumFields(*supers, Nf);
1849: for (i = 0, f = 0; i < len; ++i) {
1850: PetscInt nf, fi, ci;
1852: PetscSectionGetNumFields(s[i], &nf);
1853: for (fi = 0; fi < nf; ++fi, ++f) {
1854: const char *name = NULL;
1855: PetscInt numComp = 0;
1857: PetscSectionGetFieldName(s[i], fi, &name);
1858: PetscSectionSetFieldName(*supers, f, name);
1859: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1860: PetscSectionSetFieldComponents(*supers, f, numComp);
1861: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1862: PetscSectionGetComponentName(s[i], fi, ci, &name);
1863: PetscSectionSetComponentName(*supers, f, ci, name);
1864: }
1865: }
1866: }
1867: PetscSectionSetChart(*supers, pStart, pEnd);
1868: for (p = pStart; p < pEnd; ++p) {
1869: PetscInt dof = 0, cdof = 0;
1871: for (i = 0, f = 0; i < len; ++i) {
1872: PetscInt nf, fi, pStarti, pEndi;
1873: PetscInt fdof = 0, cfdof = 0;
1875: PetscSectionGetNumFields(s[i], &nf);
1876: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1877: if ((p < pStarti) || (p >= pEndi)) continue;
1878: for (fi = 0; fi < nf; ++fi, ++f) {
1879: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1880: PetscSectionAddFieldDof(*supers, p, f, fdof);
1881: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1882: if (cfdof) PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);
1883: dof += fdof;
1884: cdof += cfdof;
1885: }
1886: }
1887: PetscSectionSetDof(*supers, p, dof);
1888: if (cdof) PetscSectionSetConstraintDof(*supers, p, cdof);
1889: maxCdof = PetscMax(cdof, maxCdof);
1890: }
1891: PetscSectionSetUp(*supers);
1892: if (maxCdof) {
1893: PetscInt *indices;
1895: PetscMalloc1(maxCdof, &indices);
1896: for (p = pStart; p < pEnd; ++p) {
1897: PetscInt cdof;
1899: PetscSectionGetConstraintDof(*supers, p, &cdof);
1900: if (cdof) {
1901: PetscInt dof, numConst = 0, fOff = 0;
1903: for (i = 0, f = 0; i < len; ++i) {
1904: const PetscInt *oldIndices = NULL;
1905: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1907: PetscSectionGetNumFields(s[i], &nf);
1908: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1909: if ((p < pStarti) || (p >= pEndi)) continue;
1910: for (fi = 0; fi < nf; ++fi, ++f) {
1911: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1912: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1913: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1914: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
1915: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1916: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
1917: numConst += cfdof;
1918: }
1919: PetscSectionGetDof(s[i], p, &dof);
1920: fOff += dof;
1921: }
1922: PetscSectionSetConstraintIndices(*supers, p, indices);
1923: }
1924: }
1925: PetscFree(indices);
1926: }
1927: return 0;
1928: }
1930: PetscErrorCode PetscSectionCreateSubplexSection_Internal(PetscSection s, IS subpointMap, PetscBool renumberPoints, PetscSection *subs)
1931: {
1932: const PetscInt *points = NULL, *indices = NULL;
1933: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
1938: PetscSectionGetNumFields(s, &numFields);
1939: PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1940: if (numFields) PetscSectionSetNumFields(*subs, numFields);
1941: for (f = 0; f < numFields; ++f) {
1942: const char *name = NULL;
1943: PetscInt numComp = 0;
1945: PetscSectionGetFieldName(s, f, &name);
1946: PetscSectionSetFieldName(*subs, f, name);
1947: PetscSectionGetFieldComponents(s, f, &numComp);
1948: PetscSectionSetFieldComponents(*subs, f, numComp);
1949: for (c = 0; c < s->numFieldComponents[f]; ++c) {
1950: PetscSectionGetComponentName(s, f, c, &name);
1951: PetscSectionSetComponentName(*subs, f, c, name);
1952: }
1953: }
1954: /* For right now, we do not try to squeeze the subchart */
1955: if (subpointMap) {
1956: ISGetSize(subpointMap, &numSubpoints);
1957: ISGetIndices(subpointMap, &points);
1958: }
1959: PetscSectionGetChart(s, &pStart, &pEnd);
1960: if (renumberPoints) {
1961: spStart = 0;
1962: spEnd = numSubpoints;
1963: } else {
1964: ISGetMinMax(subpointMap, &spStart, &spEnd);
1965: ++spEnd;
1966: }
1967: PetscSectionSetChart(*subs, spStart, spEnd);
1968: for (p = pStart; p < pEnd; ++p) {
1969: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1971: PetscFindInt(p, numSubpoints, points, &subp);
1972: if (subp < 0) continue;
1973: if (!renumberPoints) subp = p;
1974: for (f = 0; f < numFields; ++f) {
1975: PetscSectionGetFieldDof(s, p, f, &fdof);
1976: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1977: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1978: if (cfdof) PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);
1979: }
1980: PetscSectionGetDof(s, p, &dof);
1981: PetscSectionSetDof(*subs, subp, dof);
1982: PetscSectionGetConstraintDof(s, p, &cdof);
1983: if (cdof) PetscSectionSetConstraintDof(*subs, subp, cdof);
1984: }
1985: PetscSectionSetUp(*subs);
1986: /* Change offsets to original offsets */
1987: for (p = pStart; p < pEnd; ++p) {
1988: PetscInt off, foff = 0;
1990: PetscFindInt(p, numSubpoints, points, &subp);
1991: if (subp < 0) continue;
1992: if (!renumberPoints) subp = p;
1993: for (f = 0; f < numFields; ++f) {
1994: PetscSectionGetFieldOffset(s, p, f, &foff);
1995: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1996: }
1997: PetscSectionGetOffset(s, p, &off);
1998: PetscSectionSetOffset(*subs, subp, off);
1999: }
2000: /* Copy constraint indices */
2001: for (subp = spStart; subp < spEnd; ++subp) {
2002: PetscInt cdof;
2004: PetscSectionGetConstraintDof(*subs, subp, &cdof);
2005: if (cdof) {
2006: for (f = 0; f < numFields; ++f) {
2007: PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices);
2008: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2009: }
2010: PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices);
2011: PetscSectionSetConstraintIndices(*subs, subp, indices);
2012: }
2013: }
2014: if (subpointMap) ISRestoreIndices(subpointMap, &points);
2015: return 0;
2016: }
2018: /*@
2019: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2021: Collective on s
2023: Input Parameters:
2024: + s - the `PetscSection`
2025: - subpointMap - a sorted list of points in the original mesh which are in the submesh
2027: Output Parameter:
2028: . subs - the subsection
2030: Level: advanced
2032: Note:
2033: The points are renumbered from 0, and the section offsets now refer to a new, smaller vector.
2035: Developer Note:
2036: The use of the term Submesh is confusing and unneeded
2038: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2039: @*/
2040: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2041: {
2042: PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_TRUE, subs);
2043: return 0;
2044: }
2046: /*@
2047: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2049: Collective on s
2051: Input Parameters:
2052: + s - the `PetscSection`
2053: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2055: Output Parameter:
2056: . subs - the subsection
2058: Level: advanced
2060: Note:
2061: The point numbers remain the same, but the section offsets now refer to a new, smaller vector.
2063: Developer Notes:
2064: It is unclear what the difference with `PetscSectionCreateSubmeshSection()` is.
2066: The use of the term Subdomain is unneeded and confusing
2068: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2069: @*/
2070: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2071: {
2072: PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_FALSE, subs);
2073: return 0;
2074: }
2076: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2077: {
2078: PetscInt p;
2079: PetscMPIInt rank;
2081: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2082: PetscViewerASCIIPushSynchronized(viewer);
2083: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2084: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2085: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2086: PetscInt b;
2088: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2089: if (s->bcIndices) {
2090: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]);
2091: }
2092: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2093: } else {
2094: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2095: }
2096: }
2097: PetscViewerFlush(viewer);
2098: PetscViewerASCIIPopSynchronized(viewer);
2099: if (s->sym) {
2100: PetscViewerASCIIPushTab(viewer);
2101: PetscSectionSymView(s->sym, viewer);
2102: PetscViewerASCIIPopTab(viewer);
2103: }
2104: return 0;
2105: }
2107: /*@C
2108: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2110: Collective
2112: Input Parameters:
2113: + A - the PetscSection object to view
2114: . obj - Optional object that provides the options prefix used for the options
2115: - name - command line option
2117: Level: intermediate
2119: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2120: @*/
2121: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2122: {
2124: PetscObjectViewFromOptions((PetscObject)A, obj, name);
2125: return 0;
2126: }
2128: /*@C
2129: PetscSectionView - Views a `PetscSection`
2131: Collective
2133: Input Parameters:
2134: + s - the `PetscSection` object to view
2135: - v - the viewer
2137: Level: beginner
2139: Note:
2140: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2141: distribution independent data, such as dofs, offsets, constraint dofs,
2142: and constraint indices. Points that have negative dofs, for instance,
2143: are not saved as they represent points owned by other processes.
2144: Point numbering and rank assignment is currently not stored.
2145: The saved section can be loaded with `PetscSectionLoad()`.
2147: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2148: @*/
2149: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2150: {
2151: PetscBool isascii, ishdf5;
2152: PetscInt f;
2155: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);
2157: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
2158: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2159: if (isascii) {
2160: PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer);
2161: if (s->numFields) {
2162: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields);
2163: for (f = 0; f < s->numFields; ++f) {
2164: PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
2165: PetscSectionView_ASCII(s->field[f], viewer);
2166: }
2167: } else {
2168: PetscSectionView_ASCII(s, viewer);
2169: }
2170: } else if (ishdf5) {
2171: #if PetscDefined(HAVE_HDF5)
2172: PetscSectionView_HDF5_Internal(s, viewer);
2173: #else
2174: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2175: #endif
2176: }
2177: return 0;
2178: }
2180: /*@C
2181: PetscSectionLoad - Loads a `PetscSection`
2183: Collective
2185: Input Parameters:
2186: + s - the `PetscSection` object to load
2187: - v - the viewer
2189: Level: beginner
2191: Note:
2192: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2193: a section saved with `PetscSectionView()`. The number of processes
2194: used here (N) does not need to be the same as that used when saving.
2195: After calling this function, the chart of s on rank i will be set
2196: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2197: saved section points.
2199: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2200: @*/
2201: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2202: {
2203: PetscBool ishdf5;
2207: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2208: if (ishdf5) {
2209: #if PetscDefined(HAVE_HDF5)
2210: PetscSectionLoad_HDF5_Internal(s, viewer);
2211: return 0;
2212: #else
2213: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2214: #endif
2215: } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2216: }
2218: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2219: {
2220: PetscSectionClosurePermVal clVal;
2222: if (!section->clHash) return 0;
2223: kh_foreach_value(section->clHash, clVal, {
2224: PetscFree(clVal.perm);
2225: PetscFree(clVal.invPerm);
2226: });
2227: kh_destroy(ClPerm, section->clHash);
2228: section->clHash = NULL;
2229: return 0;
2230: }
2232: /*@
2233: PetscSectionReset - Frees all section data.
2235: Not Collective
2237: Input Parameters:
2238: . s - the `PetscSection`
2240: Level: beginner
2242: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2243: @*/
2244: PetscErrorCode PetscSectionReset(PetscSection s)
2245: {
2246: PetscInt f, c;
2249: for (f = 0; f < s->numFields; ++f) {
2250: PetscSectionDestroy(&s->field[f]);
2251: PetscFree(s->fieldNames[f]);
2252: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscFree(s->compNames[f][c]);
2253: PetscFree(s->compNames[f]);
2254: }
2255: PetscFree(s->numFieldComponents);
2256: PetscFree(s->fieldNames);
2257: PetscFree(s->compNames);
2258: PetscFree(s->field);
2259: PetscSectionDestroy(&s->bc);
2260: PetscFree(s->bcIndices);
2261: PetscFree2(s->atlasDof, s->atlasOff);
2262: PetscSectionDestroy(&s->clSection);
2263: ISDestroy(&s->clPoints);
2264: ISDestroy(&s->perm);
2265: PetscSectionResetClosurePermutation(s);
2266: PetscSectionSymDestroy(&s->sym);
2267: PetscSectionDestroy(&s->clSection);
2268: ISDestroy(&s->clPoints);
2269: PetscSectionInvalidateMaxDof_Internal(s);
2270: s->pStart = -1;
2271: s->pEnd = -1;
2272: s->maxDof = 0;
2273: s->setup = PETSC_FALSE;
2274: s->numFields = 0;
2275: s->clObj = NULL;
2276: return 0;
2277: }
2279: /*@
2280: PetscSectionDestroy - Frees a section object and frees its range if that exists.
2282: Not Collective
2284: Input Parameters:
2285: . s - the `PetscSection`
2287: Level: beginner
2289: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2290: @*/
2291: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2292: {
2293: if (!*s) return 0;
2295: if (--((PetscObject)(*s))->refct > 0) {
2296: *s = NULL;
2297: return 0;
2298: }
2299: PetscSectionReset(*s);
2300: PetscHeaderDestroy(s);
2301: return 0;
2302: }
2304: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2305: {
2306: const PetscInt p = point - s->pStart;
2309: *values = &baseArray[s->atlasOff[p]];
2310: return 0;
2311: }
2313: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2314: {
2315: PetscInt *array;
2316: const PetscInt p = point - s->pStart;
2317: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2318: PetscInt cDim = 0;
2321: PetscSectionGetConstraintDof(s, p, &cDim);
2322: array = &baseArray[s->atlasOff[p]];
2323: if (!cDim) {
2324: if (orientation >= 0) {
2325: const PetscInt dim = s->atlasDof[p];
2326: PetscInt i;
2328: if (mode == INSERT_VALUES) {
2329: for (i = 0; i < dim; ++i) array[i] = values[i];
2330: } else {
2331: for (i = 0; i < dim; ++i) array[i] += values[i];
2332: }
2333: } else {
2334: PetscInt offset = 0;
2335: PetscInt j = -1, field, i;
2337: for (field = 0; field < s->numFields; ++field) {
2338: const PetscInt dim = s->field[field]->atlasDof[p];
2340: for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
2341: offset += dim;
2342: }
2343: }
2344: } else {
2345: if (orientation >= 0) {
2346: const PetscInt dim = s->atlasDof[p];
2347: PetscInt cInd = 0, i;
2348: const PetscInt *cDof;
2350: PetscSectionGetConstraintIndices(s, point, &cDof);
2351: if (mode == INSERT_VALUES) {
2352: for (i = 0; i < dim; ++i) {
2353: if ((cInd < cDim) && (i == cDof[cInd])) {
2354: ++cInd;
2355: continue;
2356: }
2357: array[i] = values[i];
2358: }
2359: } else {
2360: for (i = 0; i < dim; ++i) {
2361: if ((cInd < cDim) && (i == cDof[cInd])) {
2362: ++cInd;
2363: continue;
2364: }
2365: array[i] += values[i];
2366: }
2367: }
2368: } else {
2369: const PetscInt *cDof;
2370: PetscInt offset = 0;
2371: PetscInt cOffset = 0;
2372: PetscInt j = 0, field;
2374: PetscSectionGetConstraintIndices(s, point, &cDof);
2375: for (field = 0; field < s->numFields; ++field) {
2376: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2377: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2378: const PetscInt sDim = dim - tDim;
2379: PetscInt cInd = 0, i, k;
2381: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2382: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2383: ++cInd;
2384: continue;
2385: }
2386: array[j] = values[k];
2387: }
2388: offset += dim;
2389: cOffset += dim - tDim;
2390: }
2391: }
2392: }
2393: return 0;
2394: }
2396: /*@C
2397: PetscSectionHasConstraints - Determine whether a section has constrained dofs
2399: Not Collective
2401: Input Parameter:
2402: . s - The `PetscSection`
2404: Output Parameter:
2405: . hasConstraints - flag indicating that the section has constrained dofs
2407: Level: intermediate
2409: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2410: @*/
2411: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2412: {
2415: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2416: return 0;
2417: }
2419: /*@C
2420: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2422: Not Collective
2424: Input Parameters:
2425: + s - The `PetscSection`
2426: - point - The point
2428: Output Parameter:
2429: . indices - The constrained dofs
2431: Level: intermediate
2433: Fortran Note:
2434: Call `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`
2436: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2437: @*/
2438: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2439: {
2441: if (s->bc) {
2442: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2443: } else *indices = NULL;
2444: return 0;
2445: }
2447: /*@C
2448: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2450: Not Collective
2452: Input Parameters:
2453: + s - The `PetscSection`
2454: . point - The point
2455: - indices - The constrained dofs
2457: Level: intermediate
2459: Fortran Note:
2460: Use `PetscSectionSetConstraintIndicesF90()`
2462: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2463: @*/
2464: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2465: {
2467: if (s->bc) {
2468: const PetscInt dof = s->atlasDof[point];
2469: const PetscInt cdof = s->bc->atlasDof[point];
2470: PetscInt d;
2473: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2474: }
2475: return 0;
2476: }
2478: /*@C
2479: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2481: Not Collective
2483: Input Parameters:
2484: + s - The `PetscSection`
2485: . field - The field number
2486: - point - The point
2488: Output Parameter:
2489: . indices - The constrained dofs sorted in ascending order
2491: Level: intermediate
2493: Note:
2494: 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()`.
2496: Fortran Note:
2497: Call `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`
2499: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2500: @*/
2501: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2502: {
2505: PetscSectionCheckValidField(field, s->numFields);
2506: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2507: return 0;
2508: }
2510: /*@C
2511: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2513: Not Collective
2515: Input Parameters:
2516: + s - The `PetscSection`
2517: . point - The point
2518: . field - The field number
2519: - indices - The constrained dofs
2521: Level: intermediate
2523: Fortran Note:
2524: Use `PetscSectionSetFieldConstraintIndicesF90()`
2526: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2527: @*/
2528: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2529: {
2531: if (PetscDefined(USE_DEBUG)) {
2532: PetscInt nfdof;
2534: PetscSectionGetFieldConstraintDof(s, point, field, &nfdof);
2536: }
2537: PetscSectionCheckValidField(field, s->numFields);
2538: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2539: return 0;
2540: }
2542: /*@
2543: PetscSectionPermute - Reorder the section according to the input point permutation
2545: Collective
2547: Input Parameters:
2548: + section - The `PetscSection` object
2549: - perm - The point permutation, old point p becomes new point perm[p]
2551: Output Parameter:
2552: . sectionNew - The permuted `PetscSection`
2554: Level: intermediate
2556: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`
2557: @*/
2558: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2559: {
2560: PetscSection s = section, sNew;
2561: const PetscInt *perm;
2562: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2567: PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew);
2568: PetscSectionGetNumFields(s, &numFields);
2569: if (numFields) PetscSectionSetNumFields(sNew, numFields);
2570: for (f = 0; f < numFields; ++f) {
2571: const char *name;
2572: PetscInt numComp;
2574: PetscSectionGetFieldName(s, f, &name);
2575: PetscSectionSetFieldName(sNew, f, name);
2576: PetscSectionGetFieldComponents(s, f, &numComp);
2577: PetscSectionSetFieldComponents(sNew, f, numComp);
2578: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2579: PetscSectionGetComponentName(s, f, c, &name);
2580: PetscSectionSetComponentName(sNew, f, c, name);
2581: }
2582: }
2583: ISGetLocalSize(permutation, &numPoints);
2584: ISGetIndices(permutation, &perm);
2585: PetscSectionGetChart(s, &pStart, &pEnd);
2586: PetscSectionSetChart(sNew, pStart, pEnd);
2588: for (p = pStart; p < pEnd; ++p) {
2589: PetscInt dof, cdof;
2591: PetscSectionGetDof(s, p, &dof);
2592: PetscSectionSetDof(sNew, perm[p], dof);
2593: PetscSectionGetConstraintDof(s, p, &cdof);
2594: if (cdof) PetscSectionSetConstraintDof(sNew, perm[p], cdof);
2595: for (f = 0; f < numFields; ++f) {
2596: PetscSectionGetFieldDof(s, p, f, &dof);
2597: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2598: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2599: if (cdof) PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);
2600: }
2601: }
2602: PetscSectionSetUp(sNew);
2603: for (p = pStart; p < pEnd; ++p) {
2604: const PetscInt *cind;
2605: PetscInt cdof;
2607: PetscSectionGetConstraintDof(s, p, &cdof);
2608: if (cdof) {
2609: PetscSectionGetConstraintIndices(s, p, &cind);
2610: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2611: }
2612: for (f = 0; f < numFields; ++f) {
2613: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2614: if (cdof) {
2615: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2616: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2617: }
2618: }
2619: }
2620: ISRestoreIndices(permutation, &perm);
2621: *sectionNew = sNew;
2622: return 0;
2623: }
2625: /*@
2626: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2628: Collective
2630: Input Parameters:
2631: + section - The `PetscSection`
2632: . obj - A `PetscObject` which serves as the key for this index
2633: . clSection - `PetscSection` giving the size of the closure of each point
2634: - clPoints - `IS` giving the points in each closure
2636: Level: advanced
2638: Note:
2639: We compress out closure points with no dofs in this section
2641: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
2642: @*/
2643: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2644: {
2648: if (section->clObj != obj) PetscSectionResetClosurePermutation(section);
2649: section->clObj = obj;
2650: PetscObjectReference((PetscObject)clSection);
2651: PetscObjectReference((PetscObject)clPoints);
2652: PetscSectionDestroy(§ion->clSection);
2653: ISDestroy(§ion->clPoints);
2654: section->clSection = clSection;
2655: section->clPoints = clPoints;
2656: return 0;
2657: }
2659: /*@
2660: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
2662: Collective
2664: Input Parameters:
2665: + section - The `PetscSection`
2666: - obj - A `PetscObject` which serves as the key for this index
2668: Output Parameters:
2669: + clSection - `PetscSection` giving the size of the closure of each point
2670: - clPoints - `IS` giving the points in each closure
2672: Level: advanced
2674: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2675: @*/
2676: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2677: {
2678: if (section->clObj == obj) {
2679: if (clSection) *clSection = section->clSection;
2680: if (clPoints) *clPoints = section->clPoints;
2681: } else {
2682: if (clSection) *clSection = NULL;
2683: if (clPoints) *clPoints = NULL;
2684: }
2685: return 0;
2686: }
2688: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2689: {
2690: PetscInt i;
2691: khiter_t iter;
2692: int new_entry;
2693: PetscSectionClosurePermKey key = {depth, clSize};
2694: PetscSectionClosurePermVal *val;
2696: if (section->clObj != obj) {
2697: PetscSectionDestroy(§ion->clSection);
2698: ISDestroy(§ion->clPoints);
2699: }
2700: section->clObj = obj;
2701: if (!section->clHash) PetscClPermCreate(§ion->clHash);
2702: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2703: val = &kh_val(section->clHash, iter);
2704: if (!new_entry) {
2705: PetscFree(val->perm);
2706: PetscFree(val->invPerm);
2707: }
2708: if (mode == PETSC_COPY_VALUES) {
2709: PetscMalloc1(clSize, &val->perm);
2710: PetscArraycpy(val->perm, clPerm, clSize);
2711: } else if (mode == PETSC_OWN_POINTER) {
2712: val->perm = clPerm;
2713: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2714: PetscMalloc1(clSize, &val->invPerm);
2715: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2716: return 0;
2717: }
2719: /*@
2720: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2722: Not Collective
2724: Input Parameters:
2725: + section - The `PetscSection`
2726: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
2727: . depth - Depth of points on which to apply the given permutation
2728: - perm - Permutation of the cell dof closure
2730: Level: intermediate
2732: Notes:
2733: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2734: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2735: each topology and degree.
2737: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2738: exotic/enriched spaces on mixed topology meshes.
2740: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
2741: @*/
2742: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2743: {
2744: const PetscInt *clPerm = NULL;
2745: PetscInt clSize = 0;
2747: if (perm) {
2748: ISGetLocalSize(perm, &clSize);
2749: ISGetIndices(perm, &clPerm);
2750: }
2751: PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm);
2752: if (perm) ISRestoreIndices(perm, &clPerm);
2753: return 0;
2754: }
2756: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2757: {
2758: if (section->clObj == obj) {
2759: PetscSectionClosurePermKey k = {depth, size};
2760: PetscSectionClosurePermVal v;
2761: PetscClPermGet(section->clHash, k, &v);
2762: if (perm) *perm = v.perm;
2763: } else {
2764: if (perm) *perm = NULL;
2765: }
2766: return 0;
2767: }
2769: /*@
2770: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2772: Not Collective
2774: Input Parameters:
2775: + section - The `PetscSection`
2776: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
2777: . depth - Depth stratum on which to obtain closure permutation
2778: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2780: Output Parameter:
2781: . perm - The dof closure permutation
2783: Level: intermediate
2785: Note:
2786: The user must destroy the `IS` that is returned.
2788: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2789: @*/
2790: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2791: {
2792: const PetscInt *clPerm;
2794: PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2795: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2796: return 0;
2797: }
2799: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2800: {
2801: if (section->clObj == obj && section->clHash) {
2802: PetscSectionClosurePermKey k = {depth, size};
2803: PetscSectionClosurePermVal v;
2804: PetscClPermGet(section->clHash, k, &v);
2805: if (perm) *perm = v.invPerm;
2806: } else {
2807: if (perm) *perm = NULL;
2808: }
2809: return 0;
2810: }
2812: /*@
2813: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2815: Not Collective
2817: Input Parameters:
2818: + section - The `PetscSection`
2819: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
2820: . depth - Depth stratum on which to obtain closure permutation
2821: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2823: Output Parameters:
2824: . perm - The dof closure permutation
2826: Level: intermediate
2828: Note:
2829: The user must destroy the `IS` that is returned.
2831: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2832: @*/
2833: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2834: {
2835: const PetscInt *clPerm;
2837: PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2838: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2839: return 0;
2840: }
2842: /*@
2843: PetscSectionGetField - Get the subsection associated with a single field
2845: Input Parameters:
2846: + s - The `PetscSection`
2847: - field - The field number
2849: Output Parameter:
2850: . subs - The subsection for the given field
2852: Level: intermediate
2854: Note:
2855: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
2857: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
2858: @*/
2859: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2860: {
2863: PetscSectionCheckValidField(field, s->numFields);
2864: *subs = s->field[field];
2865: return 0;
2866: }
2868: PetscClassId PETSC_SECTION_SYM_CLASSID;
2869: PetscFunctionList PetscSectionSymList = NULL;
2871: /*@
2872: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
2874: Collective
2876: Input Parameter:
2877: . comm - the MPI communicator
2879: Output Parameter:
2880: . sym - pointer to the new set of symmetries
2882: Level: developer
2884: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
2885: @*/
2886: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2887: {
2889: ISInitializePackage();
2890: PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView);
2891: return 0;
2892: }
2894: /*@C
2895: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
2897: Collective
2899: Input Parameters:
2900: + sym - The section symmetry object
2901: - method - The name of the section symmetry type
2903: Level: developer
2905: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
2906: @*/
2907: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2908: {
2909: PetscErrorCode (*r)(PetscSectionSym);
2910: PetscBool match;
2913: PetscObjectTypeCompare((PetscObject)sym, method, &match);
2914: if (match) return 0;
2916: PetscFunctionListFind(PetscSectionSymList, method, &r);
2918: PetscTryTypeMethod(sym, destroy);
2919: sym->ops->destroy = NULL;
2921: (*r)(sym);
2922: PetscObjectChangeTypeName((PetscObject)sym, method);
2923: return 0;
2924: }
2926: /*@C
2927: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
2929: Not Collective
2931: Input Parameter:
2932: . sym - The section symmetry
2934: Output Parameter:
2935: . type - The index set type name
2937: Level: developer
2939: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
2940: @*/
2941: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2942: {
2945: *type = ((PetscObject)sym)->type_name;
2946: return 0;
2947: }
2949: /*@C
2950: PetscSectionSymRegister - Registers a new section symmetry implementation
2952: Not Collective
2954: Input Parameters:
2955: + name - The name of a new user-defined creation routine
2956: - create_func - The creation routine itself
2958: Level: developer
2960: Notes:
2961: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
2963: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
2964: @*/
2965: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2966: {
2967: ISInitializePackage();
2968: PetscFunctionListAdd(&PetscSectionSymList, sname, function);
2969: return 0;
2970: }
2972: /*@
2973: PetscSectionSymDestroy - Destroys a section symmetry.
2975: Collective
2977: Input Parameters:
2978: . sym - the section symmetry
2980: Level: developer
2982: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSymDestroy()`
2983: @*/
2984: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2985: {
2986: SymWorkLink link, next;
2988: if (!*sym) return 0;
2990: if (--((PetscObject)(*sym))->refct > 0) {
2991: *sym = NULL;
2992: return 0;
2993: }
2994: if ((*sym)->ops->destroy) (*(*sym)->ops->destroy)(*sym);
2996: for (link = (*sym)->workin; link; link = next) {
2997: PetscInt **perms = (PetscInt **)link->perms;
2998: PetscScalar **rots = (PetscScalar **)link->rots;
2999: PetscFree2(perms, rots);
3000: next = link->next;
3001: PetscFree(link);
3002: }
3003: (*sym)->workin = NULL;
3004: PetscHeaderDestroy(sym);
3005: return 0;
3006: }
3008: /*@C
3009: PetscSectionSymView - Displays a section symmetry
3011: Collective
3013: Input Parameters:
3014: + sym - the index set
3015: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3017: Level: developer
3019: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3020: @*/
3021: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3022: {
3024: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer);
3027: PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer);
3028: PetscTryTypeMethod(sym, view, viewer);
3029: return 0;
3030: }
3032: /*@
3033: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3035: Collective
3037: Input Parameters:
3038: + section - the section describing data layout
3039: - sym - the symmetry describing the affect of orientation on the access of the data
3041: Level: developer
3043: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3044: @*/
3045: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3046: {
3048: PetscSectionSymDestroy(&(section->sym));
3049: if (sym) {
3052: PetscObjectReference((PetscObject)sym);
3053: }
3054: section->sym = sym;
3055: return 0;
3056: }
3058: /*@
3059: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3061: Not Collective
3063: Input Parameters:
3064: . section - the section describing data layout
3066: Output Parameters:
3067: . sym - the symmetry describing the affect of orientation on the access of the data
3069: Level: developer
3071: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3072: @*/
3073: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3074: {
3076: *sym = section->sym;
3077: return 0;
3078: }
3080: /*@
3081: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3083: Collective
3085: Input Parameters:
3086: + section - the section describing data layout
3087: . field - the field number
3088: - sym - the symmetry describing the affect of orientation on the access of the data
3090: Level: developer
3092: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3093: @*/
3094: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3095: {
3097: PetscSectionCheckValidField(field, section->numFields);
3098: PetscSectionSetSym(section->field[field], sym);
3099: return 0;
3100: }
3102: /*@
3103: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3105: Collective
3107: Input Parameters:
3108: + section - the section describing data layout
3109: - field - the field number
3111: Output Parameters:
3112: . sym - the symmetry describing the affect of orientation on the access of the data
3114: Level: developer
3116: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3117: @*/
3118: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3119: {
3121: PetscSectionCheckValidField(field, section->numFields);
3122: *sym = section->field[field]->sym;
3123: return 0;
3124: }
3126: /*@C
3127: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3129: Not Collective
3131: Input Parameters:
3132: + section - the section
3133: . numPoints - the number of points
3134: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3135: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3136: context, see DMPlexGetConeOrientation()).
3138: Output Parameters:
3139: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3140: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3141: identity).
3143: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3144: .vb
3145: const PetscInt **perms;
3146: const PetscScalar **rots;
3147: PetscInt lOffset;
3149: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3150: for (i = 0, lOffset = 0; i < numPoints; i++) {
3151: PetscInt point = points[2*i], dof, sOffset;
3152: const PetscInt *perm = perms ? perms[i] : NULL;
3153: const PetscScalar *rot = rots ? rots[i] : NULL;
3155: PetscSectionGetDof(section,point,&dof);
3156: PetscSectionGetOffset(section,point,&sOffset);
3158: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3159: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3160: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3161: lOffset += dof;
3162: }
3163: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3164: .ve
3166: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3167: .vb
3168: const PetscInt **perms;
3169: const PetscScalar **rots;
3170: PetscInt lOffset;
3172: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3173: for (i = 0, lOffset = 0; i < numPoints; i++) {
3174: PetscInt point = points[2*i], dof, sOffset;
3175: const PetscInt *perm = perms ? perms[i] : NULL;
3176: const PetscScalar *rot = rots ? rots[i] : NULL;
3178: PetscSectionGetDof(section,point,&dof);
3179: PetscSectionGetOffset(section,point,&sOff);
3181: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3182: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3183: offset += dof;
3184: }
3185: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3186: .ve
3188: Level: developer
3190: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3191: @*/
3192: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3193: {
3194: PetscSectionSym sym;
3198: if (perms) *perms = NULL;
3199: if (rots) *rots = NULL;
3200: sym = section->sym;
3201: if (sym && (perms || rots)) {
3202: SymWorkLink link;
3204: if (sym->workin) {
3205: link = sym->workin;
3206: sym->workin = sym->workin->next;
3207: } else {
3208: PetscNew(&link);
3209: }
3210: if (numPoints > link->numPoints) {
3211: PetscInt **perms = (PetscInt **)link->perms;
3212: PetscScalar **rots = (PetscScalar **)link->rots;
3213: PetscFree2(perms, rots);
3214: PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots);
3215: link->numPoints = numPoints;
3216: }
3217: link->next = sym->workout;
3218: sym->workout = link;
3219: PetscArrayzero((PetscInt **)link->perms, numPoints);
3220: PetscArrayzero((PetscInt **)link->rots, numPoints);
3221: (*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots);
3222: if (perms) *perms = link->perms;
3223: if (rots) *rots = link->rots;
3224: }
3225: return 0;
3226: }
3228: /*@C
3229: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3231: Not Collective
3233: Input Parameters:
3234: + section - the section
3235: . numPoints - the number of points
3236: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3237: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3238: context, see `DMPlexGetConeOrientation()`).
3240: Output Parameters:
3241: + perms - The permutations for the given orientations: set to NULL at conclusion
3242: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3244: Level: developer
3246: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3247: @*/
3248: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3249: {
3250: PetscSectionSym sym;
3253: sym = section->sym;
3254: if (sym && (perms || rots)) {
3255: SymWorkLink *p, link;
3257: for (p = &sym->workout; (link = *p); p = &link->next) {
3258: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3259: *p = link->next;
3260: link->next = sym->workin;
3261: sym->workin = link;
3262: if (perms) *perms = NULL;
3263: if (rots) *rots = NULL;
3264: return 0;
3265: }
3266: }
3267: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3268: }
3269: return 0;
3270: }
3272: /*@C
3273: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3275: Not Collective
3277: Input Parameters:
3278: + section - the section
3279: . field - the field of the section
3280: . numPoints - the number of points
3281: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3282: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3283: context, see `DMPlexGetConeOrientation()`).
3285: Output Parameters:
3286: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3287: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3288: identity).
3290: Level: developer
3292: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3293: @*/
3294: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3295: {
3298: PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots);
3299: return 0;
3300: }
3302: /*@C
3303: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3305: Not Collective
3307: Input Parameters:
3308: + section - the section
3309: . field - the field number
3310: . numPoints - the number of points
3311: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3312: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3313: context, see `DMPlexGetConeOrientation()`).
3315: Output Parameters:
3316: + perms - The permutations for the given orientations: set to NULL at conclusion
3317: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3319: Level: developer
3321: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3322: @*/
3323: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3324: {
3327: PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots);
3328: return 0;
3329: }
3331: /*@
3332: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3334: Not Collective
3336: Input Parameter:
3337: . sym - the `PetscSectionSym`
3339: Output Parameter:
3340: . nsym - the equivalent symmetries
3342: Level: developer
3344: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3345: @*/
3346: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3347: {
3350: PetscTryTypeMethod(sym, copy, nsym);
3351: return 0;
3352: }
3354: /*@
3355: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3357: Collective
3359: Input Parameters:
3360: + sym - the `PetscSectionSym`
3361: - migrationSF - the distribution map from roots to leaves
3363: Output Parameters:
3364: . dsym - the redistributed symmetries
3366: Level: developer
3368: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3369: @*/
3370: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3371: {
3375: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3376: return 0;
3377: }
3379: /*@
3380: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3382: Not Collective
3384: Input Parameter:
3385: . s - the global `PetscSection`
3387: Output Parameters:
3388: . flg - the flag
3390: Level: developer
3392: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3393: @*/
3394: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3395: {
3397: *flg = s->useFieldOff;
3398: return 0;
3399: }
3401: /*@
3402: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3404: Not Collective
3406: Input Parameters:
3407: + s - the global `PetscSection`
3408: - flg - the flag
3410: Level: developer
3412: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3413: @*/
3414: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3415: {
3417: s->useFieldOff = flg;
3418: return 0;
3419: }
3421: #define PetscSectionExpandPoints_Loop(TYPE) \
3422: { \
3423: PetscInt i, n, o0, o1, size; \
3424: TYPE *a0 = (TYPE *)origArray, *a1; \
3425: PetscSectionGetStorageSize(s, &size); \
3426: PetscMalloc1(size, &a1); \
3427: for (i = 0; i < npoints; i++) { \
3428: PetscSectionGetOffset(origSection, points_[i], &o0); \
3429: PetscSectionGetOffset(s, i, &o1); \
3430: PetscSectionGetDof(s, i, &n); \
3431: PetscMemcpy(&a1[o1], &a0[o0], n *unitsize); \
3432: } \
3433: *newArray = (void *)a1; \
3434: }
3436: /*@
3437: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3439: Not Collective
3441: Input Parameters:
3442: + origSection - the `PetscSection` describing the layout of the array
3443: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
3444: . origArray - the array; its size must be equal to the storage size of origSection
3445: - points - `IS` with points to extract; its indices must lie in the chart of origSection
3447: Output Parameters:
3448: + newSection - the new `PetscSection` desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3449: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3451: Level: developer
3453: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3454: @*/
3455: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3456: {
3457: PetscSection s;
3458: const PetscInt *points_;
3459: PetscInt i, n, npoints, pStart, pEnd;
3460: PetscMPIInt unitsize;
3467: MPI_Type_size(dataType, &unitsize);
3468: ISGetLocalSize(points, &npoints);
3469: ISGetIndices(points, &points_);
3470: PetscSectionGetChart(origSection, &pStart, &pEnd);
3471: PetscSectionCreate(PETSC_COMM_SELF, &s);
3472: PetscSectionSetChart(s, 0, npoints);
3473: for (i = 0; i < npoints; i++) {
3475: PetscSectionGetDof(origSection, points_[i], &n);
3476: PetscSectionSetDof(s, i, n);
3477: }
3478: PetscSectionSetUp(s);
3479: if (newArray) {
3480: if (dataType == MPIU_INT) {
3481: PetscSectionExpandPoints_Loop(PetscInt);
3482: } else if (dataType == MPIU_SCALAR) {
3483: PetscSectionExpandPoints_Loop(PetscScalar);
3484: } else if (dataType == MPIU_REAL) {
3485: PetscSectionExpandPoints_Loop(PetscReal);
3486: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3487: }
3488: if (newSection) {
3489: *newSection = s;
3490: } else {
3491: PetscSectionDestroy(&s);
3492: }
3493: ISRestoreIndices(points, &points_);
3494: return 0;
3495: }