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. The indices returned by the `PetscSection` are appropriate for the kind of `Vec` it is associated with. For example, if the vector being indexed is a local vector, we call the section a local section. If the section indexes a global vector, we call it a global section. For parallel vectors, like global vectors, we use negative indices to indicate dofs owned by other processes.
35: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionDestroy()`, `PetscSectionCreateGlobalSection()`
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: Output 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: Note:
776: In a global section, this size will be negative for points not owned by this process.
778: Level: intermediate
780: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
781: @*/
782: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
783: {
787: 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);
788: *numDof = s->atlasDof[point - s->pStart];
789: return 0;
790: }
792: /*@
793: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
795: Not Collective
797: Input Parameters:
798: + s - the `PetscSection`
799: . point - the point
800: - numDof - the number of dof
802: Level: intermediate
804: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
805: @*/
806: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
807: {
809: 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);
810: s->atlasDof[point - s->pStart] = numDof;
811: PetscSectionInvalidateMaxDof_Internal(s);
812: return 0;
813: }
815: /*@
816: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
818: Not Collective
820: Input Parameters:
821: + s - the `PetscSection`
822: . point - the point
823: - numDof - the number of additional dof
825: Level: intermediate
827: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
828: @*/
829: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
830: {
833: 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);
834: s->atlasDof[point - s->pStart] += numDof;
835: PetscSectionInvalidateMaxDof_Internal(s);
836: return 0;
837: }
839: /*@
840: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
842: Not Collective
844: Input Parameters:
845: + s - the `PetscSection`
846: . point - the point
847: - field - the field
849: Output Parameter:
850: . numDof - the number of dof
852: Level: intermediate
854: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
855: @*/
856: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
857: {
860: PetscSectionCheckValidField(field, s->numFields);
861: PetscSectionGetDof(s->field[field], point, numDof);
862: return 0;
863: }
865: /*@
866: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
868: Not Collective
870: Input Parameters:
871: + s - the `PetscSection`
872: . point - the point
873: . field - the field
874: - numDof - the number of dof
876: Level: intermediate
878: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
879: @*/
880: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
881: {
883: PetscSectionCheckValidField(field, s->numFields);
884: PetscSectionSetDof(s->field[field], point, numDof);
885: return 0;
886: }
888: /*@
889: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
891: Not Collective
893: Input Parameters:
894: + s - the `PetscSection`
895: . point - the point
896: . field - the field
897: - numDof - the number of dof
899: Level: intermediate
901: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
902: @*/
903: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
904: {
906: PetscSectionCheckValidField(field, s->numFields);
907: PetscSectionAddDof(s->field[field], point, numDof);
908: return 0;
909: }
911: /*@
912: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
914: Not Collective
916: Input Parameters:
917: + s - the `PetscSection`
918: - point - the point
920: Output Parameter:
921: . numDof - the number of dof which are fixed by constraints
923: Level: intermediate
925: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
926: @*/
927: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
928: {
931: if (s->bc) {
932: PetscSectionGetDof(s->bc, point, numDof);
933: } else *numDof = 0;
934: return 0;
935: }
937: /*@
938: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
940: Not Collective
942: Input Parameters:
943: + s - the `PetscSection`
944: . point - the point
945: - numDof - the number of dof which are fixed by constraints
947: Level: intermediate
949: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
950: @*/
951: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
952: {
954: if (numDof) {
955: PetscSectionCheckConstraints_Private(s);
956: PetscSectionSetDof(s->bc, point, numDof);
957: }
958: return 0;
959: }
961: /*@
962: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
964: Not Collective
966: Input Parameters:
967: + s - the `PetscSection`
968: . point - the point
969: - numDof - the number of additional dof which are fixed by constraints
971: Level: intermediate
973: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
974: @*/
975: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
976: {
978: if (numDof) {
979: PetscSectionCheckConstraints_Private(s);
980: PetscSectionAddDof(s->bc, point, numDof);
981: }
982: return 0;
983: }
985: /*@
986: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
988: Not Collective
990: Input Parameters:
991: + s - the `PetscSection`
992: . point - the point
993: - field - the field
995: Output Parameter:
996: . numDof - the number of dof which are fixed by constraints
998: Level: intermediate
1000: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1001: @*/
1002: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1003: {
1006: PetscSectionCheckValidField(field, s->numFields);
1007: PetscSectionGetConstraintDof(s->field[field], point, numDof);
1008: return 0;
1009: }
1011: /*@
1012: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1014: Not Collective
1016: Input Parameters:
1017: + s - the `PetscSection`
1018: . point - the point
1019: . field - the field
1020: - numDof - the number of dof which are fixed by constraints
1022: Level: intermediate
1024: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1025: @*/
1026: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1027: {
1029: PetscSectionCheckValidField(field, s->numFields);
1030: PetscSectionSetConstraintDof(s->field[field], point, numDof);
1031: return 0;
1032: }
1034: /*@
1035: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1037: Not Collective
1039: Input Parameters:
1040: + s - the `PetscSection`
1041: . point - the point
1042: . field - the field
1043: - numDof - the number of additional dof which are fixed by constraints
1045: Level: intermediate
1047: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1048: @*/
1049: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1050: {
1052: PetscSectionCheckValidField(field, s->numFields);
1053: PetscSectionAddConstraintDof(s->field[field], point, numDof);
1054: return 0;
1055: }
1057: /*@
1058: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1060: Not Collective
1062: Input Parameter:
1063: . s - the `PetscSection`
1065: Level: advanced
1067: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1068: @*/
1069: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1070: {
1072: if (s->bc) {
1073: const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;
1075: PetscSectionSetUp(s->bc);
1076: PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices);
1077: }
1078: return 0;
1079: }
1081: /*@
1082: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
1084: Not Collective
1086: Input Parameter:
1087: . s - the `PetscSection`
1089: Level: intermediate
1091: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1092: @*/
1093: PetscErrorCode PetscSectionSetUp(PetscSection s)
1094: {
1095: const PetscInt *pind = NULL;
1096: PetscInt offset = 0, foff, p, f;
1099: if (s->setup) return 0;
1100: s->setup = PETSC_TRUE;
1101: /* Set offsets and field offsets for all points */
1102: /* Assume that all fields have the same chart */
1104: if (s->perm) ISGetIndices(s->perm, &pind);
1105: if (s->pointMajor) {
1106: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1107: const PetscInt q = pind ? pind[p] : p;
1109: /* Set point offset */
1110: s->atlasOff[q] = offset;
1111: offset += s->atlasDof[q];
1112: /* Set field offset */
1113: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1114: PetscSection sf = s->field[f];
1116: sf->atlasOff[q] = foff;
1117: foff += sf->atlasDof[q];
1118: }
1119: }
1120: } else {
1121: /* Set field offsets for all points */
1122: for (f = 0; f < s->numFields; ++f) {
1123: PetscSection sf = s->field[f];
1125: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1126: const PetscInt q = pind ? pind[p] : p;
1128: sf->atlasOff[q] = offset;
1129: offset += sf->atlasDof[q];
1130: }
1131: }
1132: /* Disable point offsets since these are unused */
1133: for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1134: }
1135: if (s->perm) ISRestoreIndices(s->perm, &pind);
1136: /* Setup BC sections */
1137: PetscSectionSetUpBC(s);
1138: for (f = 0; f < s->numFields; ++f) PetscSectionSetUpBC(s->field[f]);
1139: return 0;
1140: }
1142: /*@
1143: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1145: Not Collective
1147: Input Parameters:
1148: . s - the `PetscSection`
1150: Output Parameter:
1151: . maxDof - the maximum dof
1153: Level: intermediate
1155: Note:
1156: The returned number is up-to-date without need for `PetscSectionSetUp()`.
1158: Developer Note:
1159: The returned number is calculated lazily and stashed.
1161: A call to `PetscSectionInvalidateMaxDof_Internal()` invalidates the stashed value.
1163: `PetscSectionInvalidateMaxDof_Internal()` is called in `PetscSectionSetDof()`, `PetscSectionAddDof()` and `PetscSectionReset()`
1165: It should also be called every time atlasDof is modified directly.
1167: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1168: @*/
1169: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1170: {
1171: PetscInt p;
1175: if (s->maxDof == PETSC_MIN_INT) {
1176: s->maxDof = 0;
1177: for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1178: }
1179: *maxDof = s->maxDof;
1180: return 0;
1181: }
1183: /*@
1184: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
1186: Not Collective
1188: Input Parameter:
1189: . s - the `PetscSection`
1191: Output Parameter:
1192: . size - the size of an array which can hold all the dofs
1194: Level: intermediate
1196: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1197: @*/
1198: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1199: {
1200: PetscInt p, n = 0;
1204: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1205: *size = n;
1206: return 0;
1207: }
1209: /*@
1210: PetscSectionGetConstrainedStorageSize - Return the size of an array or local `Vec` capable of holding all unconstrained degrees of freedom.
1212: Not Collective
1214: Input Parameter:
1215: . s - the `PetscSection`
1217: Output Parameter:
1218: . size - the size of an array which can hold all unconstrained dofs
1220: Level: intermediate
1222: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1223: @*/
1224: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1225: {
1226: PetscInt p, n = 0;
1230: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1231: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1232: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1233: }
1234: *size = n;
1235: return 0;
1236: }
1238: /*@
1239: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1240: the local section and a `PetscSF` describing the section point overlap.
1242: Input Parameters:
1243: + s - The `PetscSection` for the local field layout
1244: . sf - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1245: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1246: - localOffsets - If `PETSC_TRUE`, use local rather than global offsets for the points
1248: Output Parameter:
1249: . gsection - The `PetscSection` for the global field layout
1251: Level: intermediate
1253: Notes:
1254: If we have a set of local sections defining the layout of a set of local vectors, and also a `PetscSF` to determine which section points are shared and the ownership, we can calculate a global section defining the parallel data layout, and the associated global vector.
1256: This gives negative sizes and offsets to points not owned by this process
1258: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1259: @*/
1260: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1261: {
1262: PetscSection gs;
1263: const PetscInt *pind = NULL;
1264: PetscInt *recv = NULL, *neg = NULL;
1265: PetscInt pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1266: PetscInt numFields, f, numComponents;
1274: PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs);
1275: PetscSectionGetNumFields(s, &numFields);
1276: if (numFields > 0) PetscSectionSetNumFields(gs, numFields);
1277: PetscSectionGetChart(s, &pStart, &pEnd);
1278: PetscSectionSetChart(gs, pStart, pEnd);
1279: gs->includesConstraints = includeConstraints;
1280: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1281: nlocal = nroots; /* The local/leaf space matches global/root space */
1282: /* Must allocate for all points visible to SF, which may be more than this section */
1283: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1284: PetscSFGetLeafRange(sf, NULL, &maxleaf);
1287: PetscMalloc2(nroots, &neg, nlocal, &recv);
1288: PetscArrayzero(neg, nroots);
1289: }
1290: /* Mark all local points with negative dof */
1291: for (p = pStart; p < pEnd; ++p) {
1292: PetscSectionGetDof(s, p, &dof);
1293: PetscSectionSetDof(gs, p, dof);
1294: PetscSectionGetConstraintDof(s, p, &cdof);
1295: if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(gs, p, cdof);
1296: if (neg) neg[p] = -(dof + 1);
1297: }
1298: PetscSectionSetUpBC(gs);
1299: 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]);
1300: if (nroots >= 0) {
1301: PetscArrayzero(recv, nlocal);
1302: PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1303: PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1304: for (p = pStart; p < pEnd; ++p) {
1305: if (recv[p] < 0) {
1306: gs->atlasDof[p - pStart] = recv[p];
1307: PetscSectionGetDof(s, p, &dof);
1309: }
1310: }
1311: }
1312: /* Calculate new sizes, get process offset, and calculate point offsets */
1313: if (s->perm) ISGetIndices(s->perm, &pind);
1314: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1315: const PetscInt q = pind ? pind[p] : p;
1317: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1318: gs->atlasOff[q] = off;
1319: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1320: }
1321: if (!localOffsets) {
1322: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
1323: globalOff -= off;
1324: }
1325: for (p = pStart, off = 0; p < pEnd; ++p) {
1326: gs->atlasOff[p - pStart] += globalOff;
1327: if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1328: }
1329: if (s->perm) ISRestoreIndices(s->perm, &pind);
1330: /* Put in negative offsets for ghost points */
1331: if (nroots >= 0) {
1332: PetscArrayzero(recv, nlocal);
1333: PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1334: PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE);
1335: for (p = pStart; p < pEnd; ++p) {
1336: if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1337: }
1338: }
1339: PetscFree2(neg, recv);
1340: /* Set field dofs/offsets/constraints */
1341: for (f = 0; f < numFields; ++f) {
1342: gs->field[f]->includesConstraints = includeConstraints;
1343: PetscSectionGetFieldComponents(s, f, &numComponents);
1344: PetscSectionSetFieldComponents(gs, f, numComponents);
1345: }
1346: for (p = pStart; p < pEnd; ++p) {
1347: PetscSectionGetOffset(gs, p, &off);
1348: for (f = 0, foff = off; f < numFields; ++f) {
1349: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1350: if (!includeConstraints && cdof > 0) PetscSectionSetFieldConstraintDof(gs, p, f, cdof);
1351: PetscSectionGetFieldDof(s, p, f, &dof);
1352: PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof);
1353: PetscSectionSetFieldOffset(gs, p, f, foff);
1354: PetscSectionGetFieldConstraintDof(gs, p, f, &cdof);
1355: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1356: }
1357: }
1358: for (f = 0; f < numFields; ++f) {
1359: PetscSection gfs = gs->field[f];
1361: PetscSectionSetUpBC(gfs);
1362: 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]);
1363: }
1364: gs->setup = PETSC_TRUE;
1365: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1366: *gsection = gs;
1367: return 0;
1368: }
1370: /*@
1371: PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the global field layout using
1372: the local section and an `PetscSF` describing the section point overlap.
1374: Input Parameters:
1375: + s - The `PetscSection` for the local field layout
1376: . sf - The `PetscSF` describing parallel layout of the section points
1377: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1378: . numExcludes - The number of exclusion ranges
1379: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1381: Output Parameter:
1382: . gsection - The `PetscSection` for the global field layout
1384: Level: advanced
1386: Note:
1387: This gives negative sizes and offsets to points not owned by this process
1389: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1390: @*/
1391: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1392: {
1393: const PetscInt *pind = NULL;
1394: PetscInt *neg = NULL, *tmpOff = NULL;
1395: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1400: PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection);
1401: PetscSectionGetChart(s, &pStart, &pEnd);
1402: PetscSectionSetChart(*gsection, pStart, pEnd);
1403: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1404: if (nroots >= 0) {
1406: PetscCalloc1(nroots, &neg);
1407: if (nroots > pEnd - pStart) {
1408: PetscCalloc1(nroots, &tmpOff);
1409: } else {
1410: tmpOff = &(*gsection)->atlasDof[-pStart];
1411: }
1412: }
1413: /* Mark ghost points with negative dof */
1414: for (p = pStart; p < pEnd; ++p) {
1415: for (e = 0; e < numExcludes; ++e) {
1416: if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1417: PetscSectionSetDof(*gsection, p, 0);
1418: break;
1419: }
1420: }
1421: if (e < numExcludes) continue;
1422: PetscSectionGetDof(s, p, &dof);
1423: PetscSectionSetDof(*gsection, p, dof);
1424: PetscSectionGetConstraintDof(s, p, &cdof);
1425: if (!includeConstraints && cdof > 0) PetscSectionSetConstraintDof(*gsection, p, cdof);
1426: if (neg) neg[p] = -(dof + 1);
1427: }
1428: PetscSectionSetUpBC(*gsection);
1429: if (nroots >= 0) {
1430: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1431: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1432: if (nroots > pEnd - pStart) {
1433: for (p = pStart; p < pEnd; ++p) {
1434: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1435: }
1436: }
1437: }
1438: /* Calculate new sizes, get process offset, and calculate point offsets */
1439: if (s->perm) ISGetIndices(s->perm, &pind);
1440: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1441: const PetscInt q = pind ? pind[p] : p;
1443: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1444: (*gsection)->atlasOff[q] = off;
1445: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1446: }
1447: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s));
1448: globalOff -= off;
1449: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1450: (*gsection)->atlasOff[p] += globalOff;
1451: if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1452: }
1453: if (s->perm) ISRestoreIndices(s->perm, &pind);
1454: /* Put in negative offsets for ghost points */
1455: if (nroots >= 0) {
1456: if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1457: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1458: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE);
1459: if (nroots > pEnd - pStart) {
1460: for (p = pStart; p < pEnd; ++p) {
1461: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1462: }
1463: }
1464: }
1465: if (nroots >= 0 && nroots > pEnd - pStart) PetscFree(tmpOff);
1466: PetscFree(neg);
1467: return 0;
1468: }
1470: /*@
1471: PetscSectionGetPointLayout - Get the `PetscLayout` associated with the section points
1473: Collective
1475: Input Parameters:
1476: + comm - The `MPI_Comm`
1477: - s - The `PetscSection`
1479: Output Parameter:
1480: . layout - The point layout for the section
1482: Level: advanced
1484: Note:
1485: This is usually called for the default global section.
1487: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1488: @*/
1489: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1490: {
1491: PetscInt pStart, pEnd, p, localSize = 0;
1493: PetscSectionGetChart(s, &pStart, &pEnd);
1494: for (p = pStart; p < pEnd; ++p) {
1495: PetscInt dof;
1497: PetscSectionGetDof(s, p, &dof);
1498: if (dof >= 0) ++localSize;
1499: }
1500: PetscLayoutCreate(comm, layout);
1501: PetscLayoutSetLocalSize(*layout, localSize);
1502: PetscLayoutSetBlockSize(*layout, 1);
1503: PetscLayoutSetUp(*layout);
1504: return 0;
1505: }
1507: /*@
1508: PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs.
1510: Collective
1512: Input Parameters:
1513: + comm - The `MPI_Comm`
1514: - s - The `PetscSection`
1516: Output Parameter:
1517: . layout - The dof layout for the section
1519: Level: advanced
1521: Note:
1522: This is usually called for the default global section.
1524: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1525: @*/
1526: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1527: {
1528: PetscInt pStart, pEnd, p, localSize = 0;
1532: PetscSectionGetChart(s, &pStart, &pEnd);
1533: for (p = pStart; p < pEnd; ++p) {
1534: PetscInt dof, cdof;
1536: PetscSectionGetDof(s, p, &dof);
1537: PetscSectionGetConstraintDof(s, p, &cdof);
1538: if (dof - cdof > 0) localSize += dof - cdof;
1539: }
1540: PetscLayoutCreate(comm, layout);
1541: PetscLayoutSetLocalSize(*layout, localSize);
1542: PetscLayoutSetBlockSize(*layout, 1);
1543: PetscLayoutSetUp(*layout);
1544: return 0;
1545: }
1547: /*@
1548: PetscSectionGetOffset - Return the offset into an array or `Vec` for the dof associated with the given point.
1550: Not Collective
1552: Input Parameters:
1553: + s - the `PetscSection`
1554: - point - the point
1556: Output Parameter:
1557: . offset - the offset
1559: Note:
1560: In a global section, this offset will be negative for points not owned by this process.
1562: Level: intermediate
1564: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`
1565: @*/
1566: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1567: {
1571: *offset = s->atlasOff[point - s->pStart];
1572: return 0;
1573: }
1575: /*@
1576: PetscSectionSetOffset - Set the offset into an array or `Vec` for the dof associated with the given point.
1578: Not Collective
1580: Input Parameters:
1581: + s - the `PetscSection`
1582: . point - the point
1583: - offset - the offset
1585: Level: intermediate
1587: Note:
1588: The user usually does not call this function, but uses `PetscSectionSetUp()`
1590: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1591: @*/
1592: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1593: {
1596: s->atlasOff[point - s->pStart] = offset;
1597: return 0;
1598: }
1600: /*@
1601: PetscSectionGetFieldOffset - Return the offset into an array or `Vec` for the field dof associated with the given point.
1603: Not Collective
1605: Input Parameters:
1606: + s - the `PetscSection`
1607: . point - the point
1608: - field - the field
1610: Output Parameter:
1611: . offset - the offset
1613: Note:
1614: In a global section, this offset will be negative for points not owned by this process.
1616: Level: intermediate
1618: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1619: @*/
1620: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1621: {
1624: PetscSectionCheckValidField(field, s->numFields);
1625: PetscSectionGetOffset(s->field[field], point, offset);
1626: return 0;
1627: }
1629: /*@
1630: PetscSectionSetFieldOffset - Set the offset into an array or `Vec` for the dof associated with the given field at a point.
1632: Not Collective
1634: Input Parameters:
1635: + s - the PetscSection
1636: . point - the point
1637: . field - the field
1638: - offset - the offset
1640: Level: developer
1642: Note:
1643: The user usually does not call this function, but uses `PetscSectionSetUp()`
1645: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1646: @*/
1647: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1648: {
1650: PetscSectionCheckValidField(field, s->numFields);
1651: PetscSectionSetOffset(s->field[field], point, offset);
1652: return 0;
1653: }
1655: /*@
1656: PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.
1658: Not Collective
1660: Input Parameters:
1661: + s - the `PetscSection`
1662: . point - the point
1663: - field - the field
1665: Output Parameter:
1666: . offset - the offset
1668: Level: advanced
1670: Note:
1671: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1672: this point, what number is the first dof with this field.
1674: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1675: @*/
1676: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1677: {
1678: PetscInt off, foff;
1682: PetscSectionCheckValidField(field, s->numFields);
1683: PetscSectionGetOffset(s, point, &off);
1684: PetscSectionGetOffset(s->field[field], point, &foff);
1685: *offset = foff - off;
1686: return 0;
1687: }
1689: /*@
1690: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1692: Not Collective
1694: Input Parameter:
1695: . s - the `PetscSection`
1697: Output Parameters:
1698: + start - the minimum offset
1699: - end - one more than the maximum offset
1701: Level: intermediate
1703: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1704: @*/
1705: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1706: {
1707: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1710: if (s->atlasOff) {
1711: os = s->atlasOff[0];
1712: oe = s->atlasOff[0];
1713: }
1714: PetscSectionGetChart(s, &pStart, &pEnd);
1715: for (p = 0; p < pEnd - pStart; ++p) {
1716: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1718: if (off >= 0) {
1719: os = PetscMin(os, off);
1720: oe = PetscMax(oe, off + dof);
1721: }
1722: }
1723: if (start) *start = os;
1724: if (end) *end = oe;
1725: return 0;
1726: }
1728: /*@
1729: PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields
1731: Collective
1733: Input Parameters:
1734: + s - the `PetscSection`
1735: . len - the number of subfields
1736: - fields - the subfield numbers
1738: Output Parameter:
1739: . subs - the subsection
1741: Level: advanced
1743: Notes:
1744: The section offsets now refer to a new, smaller vector.
1746: This will error if a fieldnumber is out of range
1748: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
1749: @*/
1750: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1751: {
1752: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
1754: if (!len) return 0;
1758: PetscSectionGetNumFields(s, &nF);
1760: PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1761: PetscSectionSetNumFields(*subs, len);
1762: for (f = 0; f < len; ++f) {
1763: const char *name = NULL;
1764: PetscInt numComp = 0;
1766: PetscSectionGetFieldName(s, fields[f], &name);
1767: PetscSectionSetFieldName(*subs, f, name);
1768: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1769: PetscSectionSetFieldComponents(*subs, f, numComp);
1770: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1771: PetscSectionGetComponentName(s, fields[f], c, &name);
1772: PetscSectionSetComponentName(*subs, f, c, name);
1773: }
1774: }
1775: PetscSectionGetChart(s, &pStart, &pEnd);
1776: PetscSectionSetChart(*subs, pStart, pEnd);
1777: for (p = pStart; p < pEnd; ++p) {
1778: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1780: for (f = 0; f < len; ++f) {
1781: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1782: PetscSectionSetFieldDof(*subs, p, f, fdof);
1783: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1784: if (cfdof) PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);
1785: dof += fdof;
1786: cdof += cfdof;
1787: }
1788: PetscSectionSetDof(*subs, p, dof);
1789: if (cdof) PetscSectionSetConstraintDof(*subs, p, cdof);
1790: maxCdof = PetscMax(cdof, maxCdof);
1791: }
1792: PetscSectionSetUp(*subs);
1793: if (maxCdof) {
1794: PetscInt *indices;
1796: PetscMalloc1(maxCdof, &indices);
1797: for (p = pStart; p < pEnd; ++p) {
1798: PetscInt cdof;
1800: PetscSectionGetConstraintDof(*subs, p, &cdof);
1801: if (cdof) {
1802: const PetscInt *oldIndices = NULL;
1803: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1805: for (f = 0; f < len; ++f) {
1806: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1807: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1808: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1809: PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices);
1810: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
1811: numConst += cfdof;
1812: fOff += fdof;
1813: }
1814: PetscSectionSetConstraintIndices(*subs, p, indices);
1815: }
1816: }
1817: PetscFree(indices);
1818: }
1819: return 0;
1820: }
1822: /*@
1823: PetscSectionCreateSupersection - Create a new, larger section composed of the input sections
1825: Collective
1827: Input Parameters:
1828: + s - the input sections
1829: - len - the number of input sections
1831: Output Parameter:
1832: . supers - the supersection
1834: Level: advanced
1836: Note:
1837: The section offsets now refer to a new, larger vector.
1839: Developer Note:
1840: Needs to explain how the sections are composed
1842: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
1843: @*/
1844: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1845: {
1846: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1848: if (!len) return 0;
1849: for (i = 0; i < len; ++i) {
1850: PetscInt nf, pStarti, pEndi;
1852: PetscSectionGetNumFields(s[i], &nf);
1853: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1854: pStart = PetscMin(pStart, pStarti);
1855: pEnd = PetscMax(pEnd, pEndi);
1856: Nf += nf;
1857: }
1858: PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers);
1859: PetscSectionSetNumFields(*supers, Nf);
1860: for (i = 0, f = 0; i < len; ++i) {
1861: PetscInt nf, fi, ci;
1863: PetscSectionGetNumFields(s[i], &nf);
1864: for (fi = 0; fi < nf; ++fi, ++f) {
1865: const char *name = NULL;
1866: PetscInt numComp = 0;
1868: PetscSectionGetFieldName(s[i], fi, &name);
1869: PetscSectionSetFieldName(*supers, f, name);
1870: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1871: PetscSectionSetFieldComponents(*supers, f, numComp);
1872: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1873: PetscSectionGetComponentName(s[i], fi, ci, &name);
1874: PetscSectionSetComponentName(*supers, f, ci, name);
1875: }
1876: }
1877: }
1878: PetscSectionSetChart(*supers, pStart, pEnd);
1879: for (p = pStart; p < pEnd; ++p) {
1880: PetscInt dof = 0, cdof = 0;
1882: for (i = 0, f = 0; i < len; ++i) {
1883: PetscInt nf, fi, pStarti, pEndi;
1884: PetscInt fdof = 0, cfdof = 0;
1886: PetscSectionGetNumFields(s[i], &nf);
1887: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1888: if ((p < pStarti) || (p >= pEndi)) continue;
1889: for (fi = 0; fi < nf; ++fi, ++f) {
1890: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1891: PetscSectionAddFieldDof(*supers, p, f, fdof);
1892: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1893: if (cfdof) PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);
1894: dof += fdof;
1895: cdof += cfdof;
1896: }
1897: }
1898: PetscSectionSetDof(*supers, p, dof);
1899: if (cdof) PetscSectionSetConstraintDof(*supers, p, cdof);
1900: maxCdof = PetscMax(cdof, maxCdof);
1901: }
1902: PetscSectionSetUp(*supers);
1903: if (maxCdof) {
1904: PetscInt *indices;
1906: PetscMalloc1(maxCdof, &indices);
1907: for (p = pStart; p < pEnd; ++p) {
1908: PetscInt cdof;
1910: PetscSectionGetConstraintDof(*supers, p, &cdof);
1911: if (cdof) {
1912: PetscInt dof, numConst = 0, fOff = 0;
1914: for (i = 0, f = 0; i < len; ++i) {
1915: const PetscInt *oldIndices = NULL;
1916: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1918: PetscSectionGetNumFields(s[i], &nf);
1919: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1920: if ((p < pStarti) || (p >= pEndi)) continue;
1921: for (fi = 0; fi < nf; ++fi, ++f) {
1922: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1923: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1924: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1925: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
1926: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1927: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
1928: numConst += cfdof;
1929: }
1930: PetscSectionGetDof(s[i], p, &dof);
1931: fOff += dof;
1932: }
1933: PetscSectionSetConstraintIndices(*supers, p, indices);
1934: }
1935: }
1936: PetscFree(indices);
1937: }
1938: return 0;
1939: }
1941: PetscErrorCode PetscSectionCreateSubplexSection_Internal(PetscSection s, IS subpointMap, PetscBool renumberPoints, PetscSection *subs)
1942: {
1943: const PetscInt *points = NULL, *indices = NULL;
1944: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
1949: PetscSectionGetNumFields(s, &numFields);
1950: PetscSectionCreate(PetscObjectComm((PetscObject)s), subs);
1951: if (numFields) PetscSectionSetNumFields(*subs, numFields);
1952: for (f = 0; f < numFields; ++f) {
1953: const char *name = NULL;
1954: PetscInt numComp = 0;
1956: PetscSectionGetFieldName(s, f, &name);
1957: PetscSectionSetFieldName(*subs, f, name);
1958: PetscSectionGetFieldComponents(s, f, &numComp);
1959: PetscSectionSetFieldComponents(*subs, f, numComp);
1960: for (c = 0; c < s->numFieldComponents[f]; ++c) {
1961: PetscSectionGetComponentName(s, f, c, &name);
1962: PetscSectionSetComponentName(*subs, f, c, name);
1963: }
1964: }
1965: /* For right now, we do not try to squeeze the subchart */
1966: if (subpointMap) {
1967: ISGetSize(subpointMap, &numSubpoints);
1968: ISGetIndices(subpointMap, &points);
1969: }
1970: PetscSectionGetChart(s, &pStart, &pEnd);
1971: if (renumberPoints) {
1972: spStart = 0;
1973: spEnd = numSubpoints;
1974: } else {
1975: ISGetMinMax(subpointMap, &spStart, &spEnd);
1976: ++spEnd;
1977: }
1978: PetscSectionSetChart(*subs, spStart, spEnd);
1979: for (p = pStart; p < pEnd; ++p) {
1980: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1982: PetscFindInt(p, numSubpoints, points, &subp);
1983: if (subp < 0) continue;
1984: if (!renumberPoints) subp = p;
1985: for (f = 0; f < numFields; ++f) {
1986: PetscSectionGetFieldDof(s, p, f, &fdof);
1987: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1988: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1989: if (cfdof) PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);
1990: }
1991: PetscSectionGetDof(s, p, &dof);
1992: PetscSectionSetDof(*subs, subp, dof);
1993: PetscSectionGetConstraintDof(s, p, &cdof);
1994: if (cdof) PetscSectionSetConstraintDof(*subs, subp, cdof);
1995: }
1996: PetscSectionSetUp(*subs);
1997: /* Change offsets to original offsets */
1998: for (p = pStart; p < pEnd; ++p) {
1999: PetscInt off, foff = 0;
2001: PetscFindInt(p, numSubpoints, points, &subp);
2002: if (subp < 0) continue;
2003: if (!renumberPoints) subp = p;
2004: for (f = 0; f < numFields; ++f) {
2005: PetscSectionGetFieldOffset(s, p, f, &foff);
2006: PetscSectionSetFieldOffset(*subs, subp, f, foff);
2007: }
2008: PetscSectionGetOffset(s, p, &off);
2009: PetscSectionSetOffset(*subs, subp, off);
2010: }
2011: /* Copy constraint indices */
2012: for (subp = spStart; subp < spEnd; ++subp) {
2013: PetscInt cdof;
2015: PetscSectionGetConstraintDof(*subs, subp, &cdof);
2016: if (cdof) {
2017: for (f = 0; f < numFields; ++f) {
2018: PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices);
2019: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2020: }
2021: PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices);
2022: PetscSectionSetConstraintIndices(*subs, subp, indices);
2023: }
2024: }
2025: if (subpointMap) ISRestoreIndices(subpointMap, &points);
2026: return 0;
2027: }
2029: /*@
2030: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2032: Collective
2034: Input Parameters:
2035: + s - the `PetscSection`
2036: - subpointMap - a sorted list of points in the original mesh which are in the submesh
2038: Output Parameter:
2039: . subs - the subsection
2041: Level: advanced
2043: Note:
2044: The points are renumbered from 0, and the section offsets now refer to a new, smaller vector.
2046: Developer Note:
2047: The use of the term Submesh is confusing and unneeded
2049: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2050: @*/
2051: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2052: {
2053: PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_TRUE, subs);
2054: return 0;
2055: }
2057: /*@
2058: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2060: Collective
2062: Input Parameters:
2063: + s - the `PetscSection`
2064: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2066: Output Parameter:
2067: . subs - the subsection
2069: Level: advanced
2071: Note:
2072: The point numbers remain the same, but the section offsets now refer to a new, smaller vector.
2074: Developer Notes:
2075: It is unclear what the difference with `PetscSectionCreateSubmeshSection()` is.
2077: The use of the term Subdomain is unneeded and confusing
2079: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2080: @*/
2081: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2082: {
2083: PetscSectionCreateSubplexSection_Internal(s, subpointMap, PETSC_FALSE, subs);
2084: return 0;
2085: }
2087: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2088: {
2089: PetscInt p;
2090: PetscMPIInt rank;
2092: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2093: PetscViewerASCIIPushSynchronized(viewer);
2094: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2095: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2096: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2097: PetscInt b;
2099: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2100: if (s->bcIndices) {
2101: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]);
2102: }
2103: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2104: } else {
2105: PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]);
2106: }
2107: }
2108: PetscViewerFlush(viewer);
2109: PetscViewerASCIIPopSynchronized(viewer);
2110: if (s->sym) {
2111: PetscViewerASCIIPushTab(viewer);
2112: PetscSectionSymView(s->sym, viewer);
2113: PetscViewerASCIIPopTab(viewer);
2114: }
2115: return 0;
2116: }
2118: /*@C
2119: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2121: Collective
2123: Input Parameters:
2124: + A - the PetscSection object to view
2125: . obj - Optional object that provides the options prefix used for the options
2126: - name - command line option
2128: Level: intermediate
2130: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2131: @*/
2132: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2133: {
2135: PetscObjectViewFromOptions((PetscObject)A, obj, name);
2136: return 0;
2137: }
2139: /*@C
2140: PetscSectionView - Views a `PetscSection`
2142: Collective
2144: Input Parameters:
2145: + s - the `PetscSection` object to view
2146: - v - the viewer
2148: Level: beginner
2150: Note:
2151: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2152: distribution independent data, such as dofs, offsets, constraint dofs,
2153: and constraint indices. Points that have negative dofs, for instance,
2154: are not saved as they represent points owned by other processes.
2155: Point numbering and rank assignment is currently not stored.
2156: The saved section can be loaded with `PetscSectionLoad()`.
2158: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2159: @*/
2160: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2161: {
2162: PetscBool isascii, ishdf5;
2163: PetscInt f;
2166: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);
2168: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
2169: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2170: if (isascii) {
2171: PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer);
2172: if (s->numFields) {
2173: PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields);
2174: for (f = 0; f < s->numFields; ++f) {
2175: PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]);
2176: PetscSectionView_ASCII(s->field[f], viewer);
2177: }
2178: } else {
2179: PetscSectionView_ASCII(s, viewer);
2180: }
2181: } else if (ishdf5) {
2182: #if PetscDefined(HAVE_HDF5)
2183: PetscSectionView_HDF5_Internal(s, viewer);
2184: #else
2185: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2186: #endif
2187: }
2188: return 0;
2189: }
2191: /*@C
2192: PetscSectionLoad - Loads a `PetscSection`
2194: Collective
2196: Input Parameters:
2197: + s - the `PetscSection` object to load
2198: - v - the viewer
2200: Level: beginner
2202: Note:
2203: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2204: a section saved with `PetscSectionView()`. The number of processes
2205: used here (N) does not need to be the same as that used when saving.
2206: After calling this function, the chart of s on rank i will be set
2207: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2208: saved section points.
2210: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2211: @*/
2212: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2213: {
2214: PetscBool ishdf5;
2218: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
2219: if (ishdf5) {
2220: #if PetscDefined(HAVE_HDF5)
2221: PetscSectionLoad_HDF5_Internal(s, viewer);
2222: return 0;
2223: #else
2224: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2225: #endif
2226: } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2227: }
2229: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2230: {
2231: PetscSectionClosurePermVal clVal;
2233: if (!section->clHash) return 0;
2234: kh_foreach_value(section->clHash, clVal, {
2235: PetscFree(clVal.perm);
2236: PetscFree(clVal.invPerm);
2237: });
2238: kh_destroy(ClPerm, section->clHash);
2239: section->clHash = NULL;
2240: return 0;
2241: }
2243: /*@
2244: PetscSectionReset - Frees all section data.
2246: Not Collective
2248: Input Parameters:
2249: . s - the `PetscSection`
2251: Level: beginner
2253: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2254: @*/
2255: PetscErrorCode PetscSectionReset(PetscSection s)
2256: {
2257: PetscInt f, c;
2260: for (f = 0; f < s->numFields; ++f) {
2261: PetscSectionDestroy(&s->field[f]);
2262: PetscFree(s->fieldNames[f]);
2263: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscFree(s->compNames[f][c]);
2264: PetscFree(s->compNames[f]);
2265: }
2266: PetscFree(s->numFieldComponents);
2267: PetscFree(s->fieldNames);
2268: PetscFree(s->compNames);
2269: PetscFree(s->field);
2270: PetscSectionDestroy(&s->bc);
2271: PetscFree(s->bcIndices);
2272: PetscFree2(s->atlasDof, s->atlasOff);
2273: PetscSectionDestroy(&s->clSection);
2274: ISDestroy(&s->clPoints);
2275: ISDestroy(&s->perm);
2276: PetscSectionResetClosurePermutation(s);
2277: PetscSectionSymDestroy(&s->sym);
2278: PetscSectionDestroy(&s->clSection);
2279: ISDestroy(&s->clPoints);
2280: PetscSectionInvalidateMaxDof_Internal(s);
2281: s->pStart = -1;
2282: s->pEnd = -1;
2283: s->maxDof = 0;
2284: s->setup = PETSC_FALSE;
2285: s->numFields = 0;
2286: s->clObj = NULL;
2287: return 0;
2288: }
2290: /*@
2291: PetscSectionDestroy - Frees a section object and frees its range if that exists.
2293: Not Collective
2295: Input Parameters:
2296: . s - the `PetscSection`
2298: Level: beginner
2300: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2301: @*/
2302: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2303: {
2304: if (!*s) return 0;
2306: if (--((PetscObject)(*s))->refct > 0) {
2307: *s = NULL;
2308: return 0;
2309: }
2310: PetscSectionReset(*s);
2311: PetscHeaderDestroy(s);
2312: return 0;
2313: }
2315: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2316: {
2317: const PetscInt p = point - s->pStart;
2320: *values = &baseArray[s->atlasOff[p]];
2321: return 0;
2322: }
2324: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2325: {
2326: PetscInt *array;
2327: const PetscInt p = point - s->pStart;
2328: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2329: PetscInt cDim = 0;
2332: PetscSectionGetConstraintDof(s, p, &cDim);
2333: array = &baseArray[s->atlasOff[p]];
2334: if (!cDim) {
2335: if (orientation >= 0) {
2336: const PetscInt dim = s->atlasDof[p];
2337: PetscInt i;
2339: if (mode == INSERT_VALUES) {
2340: for (i = 0; i < dim; ++i) array[i] = values[i];
2341: } else {
2342: for (i = 0; i < dim; ++i) array[i] += values[i];
2343: }
2344: } else {
2345: PetscInt offset = 0;
2346: PetscInt j = -1, field, i;
2348: for (field = 0; field < s->numFields; ++field) {
2349: const PetscInt dim = s->field[field]->atlasDof[p];
2351: for (i = dim - 1; i >= 0; --i) array[++j] = values[i + offset];
2352: offset += dim;
2353: }
2354: }
2355: } else {
2356: if (orientation >= 0) {
2357: const PetscInt dim = s->atlasDof[p];
2358: PetscInt cInd = 0, i;
2359: const PetscInt *cDof;
2361: PetscSectionGetConstraintIndices(s, point, &cDof);
2362: if (mode == INSERT_VALUES) {
2363: for (i = 0; i < dim; ++i) {
2364: if ((cInd < cDim) && (i == cDof[cInd])) {
2365: ++cInd;
2366: continue;
2367: }
2368: array[i] = values[i];
2369: }
2370: } else {
2371: for (i = 0; i < dim; ++i) {
2372: if ((cInd < cDim) && (i == cDof[cInd])) {
2373: ++cInd;
2374: continue;
2375: }
2376: array[i] += values[i];
2377: }
2378: }
2379: } else {
2380: const PetscInt *cDof;
2381: PetscInt offset = 0;
2382: PetscInt cOffset = 0;
2383: PetscInt j = 0, field;
2385: PetscSectionGetConstraintIndices(s, point, &cDof);
2386: for (field = 0; field < s->numFields; ++field) {
2387: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2388: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2389: const PetscInt sDim = dim - tDim;
2390: PetscInt cInd = 0, i, k;
2392: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2393: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2394: ++cInd;
2395: continue;
2396: }
2397: array[j] = values[k];
2398: }
2399: offset += dim;
2400: cOffset += dim - tDim;
2401: }
2402: }
2403: }
2404: return 0;
2405: }
2407: /*@C
2408: PetscSectionHasConstraints - Determine whether a section has constrained dofs
2410: Not Collective
2412: Input Parameter:
2413: . s - The `PetscSection`
2415: Output Parameter:
2416: . hasConstraints - flag indicating that the section has constrained dofs
2418: Level: intermediate
2420: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2421: @*/
2422: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2423: {
2426: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2427: return 0;
2428: }
2430: /*@C
2431: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2433: Not Collective
2435: Input Parameters:
2436: + s - The `PetscSection`
2437: - point - The point
2439: Output Parameter:
2440: . indices - The constrained dofs
2442: Level: intermediate
2444: Fortran Note:
2445: Call `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`
2447: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2448: @*/
2449: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2450: {
2452: if (s->bc) {
2453: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2454: } else *indices = NULL;
2455: return 0;
2456: }
2458: /*@C
2459: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2461: Not Collective
2463: Input Parameters:
2464: + s - The `PetscSection`
2465: . point - The point
2466: - indices - The constrained dofs
2468: Level: intermediate
2470: Fortran Note:
2471: Use `PetscSectionSetConstraintIndicesF90()`
2473: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2474: @*/
2475: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2476: {
2478: if (s->bc) {
2479: const PetscInt dof = s->atlasDof[point];
2480: const PetscInt cdof = s->bc->atlasDof[point];
2481: PetscInt d;
2484: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2485: }
2486: return 0;
2487: }
2489: /*@C
2490: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2492: Not Collective
2494: Input Parameters:
2495: + s - The `PetscSection`
2496: . field - The field number
2497: - point - The point
2499: Output Parameter:
2500: . indices - The constrained dofs sorted in ascending order
2502: Level: intermediate
2504: Note:
2505: 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()`.
2507: Fortran Note:
2508: Call `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`
2510: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2511: @*/
2512: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2513: {
2516: PetscSectionCheckValidField(field, s->numFields);
2517: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2518: return 0;
2519: }
2521: /*@C
2522: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2524: Not Collective
2526: Input Parameters:
2527: + s - The `PetscSection`
2528: . point - The point
2529: . field - The field number
2530: - indices - The constrained dofs
2532: Level: intermediate
2534: Fortran Note:
2535: Use `PetscSectionSetFieldConstraintIndicesF90()`
2537: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2538: @*/
2539: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2540: {
2542: if (PetscDefined(USE_DEBUG)) {
2543: PetscInt nfdof;
2545: PetscSectionGetFieldConstraintDof(s, point, field, &nfdof);
2547: }
2548: PetscSectionCheckValidField(field, s->numFields);
2549: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2550: return 0;
2551: }
2553: /*@
2554: PetscSectionPermute - Reorder the section according to the input point permutation
2556: Collective
2558: Input Parameters:
2559: + section - The `PetscSection` object
2560: - perm - The point permutation, old point p becomes new point perm[p]
2562: Output Parameter:
2563: . sectionNew - The permuted `PetscSection`
2565: Level: intermediate
2567: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`
2568: @*/
2569: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2570: {
2571: PetscSection s = section, sNew;
2572: const PetscInt *perm;
2573: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2578: PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew);
2579: PetscSectionGetNumFields(s, &numFields);
2580: if (numFields) PetscSectionSetNumFields(sNew, numFields);
2581: for (f = 0; f < numFields; ++f) {
2582: const char *name;
2583: PetscInt numComp;
2585: PetscSectionGetFieldName(s, f, &name);
2586: PetscSectionSetFieldName(sNew, f, name);
2587: PetscSectionGetFieldComponents(s, f, &numComp);
2588: PetscSectionSetFieldComponents(sNew, f, numComp);
2589: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2590: PetscSectionGetComponentName(s, f, c, &name);
2591: PetscSectionSetComponentName(sNew, f, c, name);
2592: }
2593: }
2594: ISGetLocalSize(permutation, &numPoints);
2595: ISGetIndices(permutation, &perm);
2596: PetscSectionGetChart(s, &pStart, &pEnd);
2597: PetscSectionSetChart(sNew, pStart, pEnd);
2599: for (p = pStart; p < pEnd; ++p) {
2600: PetscInt dof, cdof;
2602: PetscSectionGetDof(s, p, &dof);
2603: PetscSectionSetDof(sNew, perm[p], dof);
2604: PetscSectionGetConstraintDof(s, p, &cdof);
2605: if (cdof) PetscSectionSetConstraintDof(sNew, perm[p], cdof);
2606: for (f = 0; f < numFields; ++f) {
2607: PetscSectionGetFieldDof(s, p, f, &dof);
2608: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2609: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2610: if (cdof) PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);
2611: }
2612: }
2613: PetscSectionSetUp(sNew);
2614: for (p = pStart; p < pEnd; ++p) {
2615: const PetscInt *cind;
2616: PetscInt cdof;
2618: PetscSectionGetConstraintDof(s, p, &cdof);
2619: if (cdof) {
2620: PetscSectionGetConstraintIndices(s, p, &cind);
2621: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2622: }
2623: for (f = 0; f < numFields; ++f) {
2624: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2625: if (cdof) {
2626: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2627: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2628: }
2629: }
2630: }
2631: ISRestoreIndices(permutation, &perm);
2632: *sectionNew = sNew;
2633: return 0;
2634: }
2636: /*@
2637: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2639: Collective
2641: Input Parameters:
2642: + section - The `PetscSection`
2643: . obj - A `PetscObject` which serves as the key for this index
2644: . clSection - `PetscSection` giving the size of the closure of each point
2645: - clPoints - `IS` giving the points in each closure
2647: Level: advanced
2649: Note:
2650: We compress out closure points with no dofs in this section
2652: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
2653: @*/
2654: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2655: {
2659: if (section->clObj != obj) PetscSectionResetClosurePermutation(section);
2660: section->clObj = obj;
2661: PetscObjectReference((PetscObject)clSection);
2662: PetscObjectReference((PetscObject)clPoints);
2663: PetscSectionDestroy(§ion->clSection);
2664: ISDestroy(§ion->clPoints);
2665: section->clSection = clSection;
2666: section->clPoints = clPoints;
2667: return 0;
2668: }
2670: /*@
2671: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
2673: Collective
2675: Input Parameters:
2676: + section - The `PetscSection`
2677: - obj - A `PetscObject` which serves as the key for this index
2679: Output Parameters:
2680: + clSection - `PetscSection` giving the size of the closure of each point
2681: - clPoints - `IS` giving the points in each closure
2683: Level: advanced
2685: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2686: @*/
2687: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2688: {
2689: if (section->clObj == obj) {
2690: if (clSection) *clSection = section->clSection;
2691: if (clPoints) *clPoints = section->clPoints;
2692: } else {
2693: if (clSection) *clSection = NULL;
2694: if (clPoints) *clPoints = NULL;
2695: }
2696: return 0;
2697: }
2699: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2700: {
2701: PetscInt i;
2702: khiter_t iter;
2703: int new_entry;
2704: PetscSectionClosurePermKey key = {depth, clSize};
2705: PetscSectionClosurePermVal *val;
2707: if (section->clObj != obj) {
2708: PetscSectionDestroy(§ion->clSection);
2709: ISDestroy(§ion->clPoints);
2710: }
2711: section->clObj = obj;
2712: if (!section->clHash) PetscClPermCreate(§ion->clHash);
2713: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2714: val = &kh_val(section->clHash, iter);
2715: if (!new_entry) {
2716: PetscFree(val->perm);
2717: PetscFree(val->invPerm);
2718: }
2719: if (mode == PETSC_COPY_VALUES) {
2720: PetscMalloc1(clSize, &val->perm);
2721: PetscArraycpy(val->perm, clPerm, clSize);
2722: } else if (mode == PETSC_OWN_POINTER) {
2723: val->perm = clPerm;
2724: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2725: PetscMalloc1(clSize, &val->invPerm);
2726: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2727: return 0;
2728: }
2730: /*@
2731: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2733: Not Collective
2735: Input Parameters:
2736: + section - The `PetscSection`
2737: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
2738: . depth - Depth of points on which to apply the given permutation
2739: - perm - Permutation of the cell dof closure
2741: Level: intermediate
2743: Notes:
2744: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2745: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2746: each topology and degree.
2748: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2749: exotic/enriched spaces on mixed topology meshes.
2751: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
2752: @*/
2753: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2754: {
2755: const PetscInt *clPerm = NULL;
2756: PetscInt clSize = 0;
2758: if (perm) {
2759: ISGetLocalSize(perm, &clSize);
2760: ISGetIndices(perm, &clPerm);
2761: }
2762: PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm);
2763: if (perm) ISRestoreIndices(perm, &clPerm);
2764: return 0;
2765: }
2767: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2768: {
2769: if (section->clObj == obj) {
2770: PetscSectionClosurePermKey k = {depth, size};
2771: PetscSectionClosurePermVal v;
2772: PetscClPermGet(section->clHash, k, &v);
2773: if (perm) *perm = v.perm;
2774: } else {
2775: if (perm) *perm = NULL;
2776: }
2777: return 0;
2778: }
2780: /*@
2781: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2783: Not Collective
2785: Input Parameters:
2786: + section - The `PetscSection`
2787: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
2788: . depth - Depth stratum on which to obtain closure permutation
2789: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2791: Output Parameter:
2792: . perm - The dof closure permutation
2794: Level: intermediate
2796: Note:
2797: The user must destroy the `IS` that is returned.
2799: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2800: @*/
2801: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2802: {
2803: const PetscInt *clPerm;
2805: PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2806: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2807: return 0;
2808: }
2810: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2811: {
2812: if (section->clObj == obj && section->clHash) {
2813: PetscSectionClosurePermKey k = {depth, size};
2814: PetscSectionClosurePermVal v;
2815: PetscClPermGet(section->clHash, k, &v);
2816: if (perm) *perm = v.invPerm;
2817: } else {
2818: if (perm) *perm = NULL;
2819: }
2820: return 0;
2821: }
2823: /*@
2824: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2826: Not Collective
2828: Input Parameters:
2829: + section - The `PetscSection`
2830: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
2831: . depth - Depth stratum on which to obtain closure permutation
2832: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2834: Output Parameters:
2835: . perm - The dof closure permutation
2837: Level: intermediate
2839: Note:
2840: The user must destroy the `IS` that is returned.
2842: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
2843: @*/
2844: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2845: {
2846: const PetscInt *clPerm;
2848: PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2849: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2850: return 0;
2851: }
2853: /*@
2854: PetscSectionGetField - Get the subsection associated with a single field
2856: Input Parameters:
2857: + s - The `PetscSection`
2858: - field - The field number
2860: Output Parameter:
2861: . subs - The subsection for the given field
2863: Level: intermediate
2865: Note:
2866: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
2868: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
2869: @*/
2870: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2871: {
2874: PetscSectionCheckValidField(field, s->numFields);
2875: *subs = s->field[field];
2876: return 0;
2877: }
2879: PetscClassId PETSC_SECTION_SYM_CLASSID;
2880: PetscFunctionList PetscSectionSymList = NULL;
2882: /*@
2883: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
2885: Collective
2887: Input Parameter:
2888: . comm - the MPI communicator
2890: Output Parameter:
2891: . sym - pointer to the new set of symmetries
2893: Level: developer
2895: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
2896: @*/
2897: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2898: {
2900: ISInitializePackage();
2901: PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView);
2902: return 0;
2903: }
2905: /*@C
2906: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
2908: Collective
2910: Input Parameters:
2911: + sym - The section symmetry object
2912: - method - The name of the section symmetry type
2914: Level: developer
2916: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
2917: @*/
2918: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2919: {
2920: PetscErrorCode (*r)(PetscSectionSym);
2921: PetscBool match;
2924: PetscObjectTypeCompare((PetscObject)sym, method, &match);
2925: if (match) return 0;
2927: PetscFunctionListFind(PetscSectionSymList, method, &r);
2929: PetscTryTypeMethod(sym, destroy);
2930: sym->ops->destroy = NULL;
2932: (*r)(sym);
2933: PetscObjectChangeTypeName((PetscObject)sym, method);
2934: return 0;
2935: }
2937: /*@C
2938: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
2940: Not Collective
2942: Input Parameter:
2943: . sym - The section symmetry
2945: Output Parameter:
2946: . type - The index set type name
2948: Level: developer
2950: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
2951: @*/
2952: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2953: {
2956: *type = ((PetscObject)sym)->type_name;
2957: return 0;
2958: }
2960: /*@C
2961: PetscSectionSymRegister - Registers a new section symmetry implementation
2963: Not Collective
2965: Input Parameters:
2966: + name - The name of a new user-defined creation routine
2967: - create_func - The creation routine itself
2969: Level: developer
2971: Notes:
2972: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
2974: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
2975: @*/
2976: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2977: {
2978: ISInitializePackage();
2979: PetscFunctionListAdd(&PetscSectionSymList, sname, function);
2980: return 0;
2981: }
2983: /*@
2984: PetscSectionSymDestroy - Destroys a section symmetry.
2986: Collective
2988: Input Parameters:
2989: . sym - the section symmetry
2991: Level: developer
2993: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSymDestroy()`
2994: @*/
2995: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2996: {
2997: SymWorkLink link, next;
2999: if (!*sym) return 0;
3001: if (--((PetscObject)(*sym))->refct > 0) {
3002: *sym = NULL;
3003: return 0;
3004: }
3005: if ((*sym)->ops->destroy) (*(*sym)->ops->destroy)(*sym);
3007: for (link = (*sym)->workin; link; link = next) {
3008: PetscInt **perms = (PetscInt **)link->perms;
3009: PetscScalar **rots = (PetscScalar **)link->rots;
3010: PetscFree2(perms, rots);
3011: next = link->next;
3012: PetscFree(link);
3013: }
3014: (*sym)->workin = NULL;
3015: PetscHeaderDestroy(sym);
3016: return 0;
3017: }
3019: /*@C
3020: PetscSectionSymView - Displays a section symmetry
3022: Collective
3024: Input Parameters:
3025: + sym - the index set
3026: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3028: Level: developer
3030: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3031: @*/
3032: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3033: {
3035: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer);
3038: PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer);
3039: PetscTryTypeMethod(sym, view, viewer);
3040: return 0;
3041: }
3043: /*@
3044: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3046: Collective
3048: Input Parameters:
3049: + section - the section describing data layout
3050: - sym - the symmetry describing the affect of orientation on the access of the data
3052: Level: developer
3054: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3055: @*/
3056: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3057: {
3059: PetscSectionSymDestroy(&(section->sym));
3060: if (sym) {
3063: PetscObjectReference((PetscObject)sym);
3064: }
3065: section->sym = sym;
3066: return 0;
3067: }
3069: /*@
3070: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3072: Not Collective
3074: Input Parameters:
3075: . section - the section describing data layout
3077: Output Parameters:
3078: . sym - the symmetry describing the affect of orientation on the access of the data
3080: Level: developer
3082: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3083: @*/
3084: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3085: {
3087: *sym = section->sym;
3088: return 0;
3089: }
3091: /*@
3092: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3094: Collective
3096: Input Parameters:
3097: + section - the section describing data layout
3098: . field - the field number
3099: - sym - the symmetry describing the affect of orientation on the access of the data
3101: Level: developer
3103: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3104: @*/
3105: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3106: {
3108: PetscSectionCheckValidField(field, section->numFields);
3109: PetscSectionSetSym(section->field[field], sym);
3110: return 0;
3111: }
3113: /*@
3114: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3116: Collective
3118: Input Parameters:
3119: + section - the section describing data layout
3120: - field - the field number
3122: Output Parameters:
3123: . sym - the symmetry describing the affect of orientation on the access of the data
3125: Level: developer
3127: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3128: @*/
3129: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3130: {
3132: PetscSectionCheckValidField(field, section->numFields);
3133: *sym = section->field[field]->sym;
3134: return 0;
3135: }
3137: /*@C
3138: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3140: Not Collective
3142: Input Parameters:
3143: + section - the section
3144: . numPoints - the number of points
3145: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3146: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3147: context, see DMPlexGetConeOrientation()).
3149: Output Parameters:
3150: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3151: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3152: identity).
3154: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3155: .vb
3156: const PetscInt **perms;
3157: const PetscScalar **rots;
3158: PetscInt lOffset;
3160: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3161: for (i = 0, lOffset = 0; i < numPoints; i++) {
3162: PetscInt point = points[2*i], dof, sOffset;
3163: const PetscInt *perm = perms ? perms[i] : NULL;
3164: const PetscScalar *rot = rots ? rots[i] : NULL;
3166: PetscSectionGetDof(section,point,&dof);
3167: PetscSectionGetOffset(section,point,&sOffset);
3169: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3170: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3171: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3172: lOffset += dof;
3173: }
3174: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3175: .ve
3177: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3178: .vb
3179: const PetscInt **perms;
3180: const PetscScalar **rots;
3181: PetscInt lOffset;
3183: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3184: for (i = 0, lOffset = 0; i < numPoints; i++) {
3185: PetscInt point = points[2*i], dof, sOffset;
3186: const PetscInt *perm = perms ? perms[i] : NULL;
3187: const PetscScalar *rot = rots ? rots[i] : NULL;
3189: PetscSectionGetDof(section,point,&dof);
3190: PetscSectionGetOffset(section,point,&sOff);
3192: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3193: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3194: offset += dof;
3195: }
3196: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3197: .ve
3199: Level: developer
3201: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3202: @*/
3203: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3204: {
3205: PetscSectionSym sym;
3209: if (perms) *perms = NULL;
3210: if (rots) *rots = NULL;
3211: sym = section->sym;
3212: if (sym && (perms || rots)) {
3213: SymWorkLink link;
3215: if (sym->workin) {
3216: link = sym->workin;
3217: sym->workin = sym->workin->next;
3218: } else {
3219: PetscNew(&link);
3220: }
3221: if (numPoints > link->numPoints) {
3222: PetscInt **perms = (PetscInt **)link->perms;
3223: PetscScalar **rots = (PetscScalar **)link->rots;
3224: PetscFree2(perms, rots);
3225: PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots);
3226: link->numPoints = numPoints;
3227: }
3228: link->next = sym->workout;
3229: sym->workout = link;
3230: PetscArrayzero((PetscInt **)link->perms, numPoints);
3231: PetscArrayzero((PetscInt **)link->rots, numPoints);
3232: (*sym->ops->getpoints)(sym, section, numPoints, points, link->perms, link->rots);
3233: if (perms) *perms = link->perms;
3234: if (rots) *rots = link->rots;
3235: }
3236: return 0;
3237: }
3239: /*@C
3240: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3242: Not Collective
3244: Input Parameters:
3245: + section - the section
3246: . numPoints - the number of points
3247: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3248: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3249: context, see `DMPlexGetConeOrientation()`).
3251: Output Parameters:
3252: + perms - The permutations for the given orientations: set to NULL at conclusion
3253: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3255: Level: developer
3257: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3258: @*/
3259: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3260: {
3261: PetscSectionSym sym;
3264: sym = section->sym;
3265: if (sym && (perms || rots)) {
3266: SymWorkLink *p, link;
3268: for (p = &sym->workout; (link = *p); p = &link->next) {
3269: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3270: *p = link->next;
3271: link->next = sym->workin;
3272: sym->workin = link;
3273: if (perms) *perms = NULL;
3274: if (rots) *rots = NULL;
3275: return 0;
3276: }
3277: }
3278: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3279: }
3280: return 0;
3281: }
3283: /*@C
3284: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3286: Not Collective
3288: Input Parameters:
3289: + section - the section
3290: . field - the field of the section
3291: . numPoints - the number of points
3292: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3293: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3294: context, see `DMPlexGetConeOrientation()`).
3296: Output Parameters:
3297: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3298: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3299: identity).
3301: Level: developer
3303: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3304: @*/
3305: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3306: {
3309: PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots);
3310: return 0;
3311: }
3313: /*@C
3314: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3316: Not Collective
3318: Input Parameters:
3319: + section - the section
3320: . field - the field number
3321: . numPoints - the number of points
3322: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3323: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3324: context, see `DMPlexGetConeOrientation()`).
3326: Output Parameters:
3327: + perms - The permutations for the given orientations: set to NULL at conclusion
3328: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3330: Level: developer
3332: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3333: @*/
3334: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3335: {
3338: PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots);
3339: return 0;
3340: }
3342: /*@
3343: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3345: Not Collective
3347: Input Parameter:
3348: . sym - the `PetscSectionSym`
3350: Output Parameter:
3351: . nsym - the equivalent symmetries
3353: Level: developer
3355: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3356: @*/
3357: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3358: {
3361: PetscTryTypeMethod(sym, copy, nsym);
3362: return 0;
3363: }
3365: /*@
3366: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3368: Collective
3370: Input Parameters:
3371: + sym - the `PetscSectionSym`
3372: - migrationSF - the distribution map from roots to leaves
3374: Output Parameters:
3375: . dsym - the redistributed symmetries
3377: Level: developer
3379: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3380: @*/
3381: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3382: {
3386: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3387: return 0;
3388: }
3390: /*@
3391: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3393: Not Collective
3395: Input Parameter:
3396: . s - the global `PetscSection`
3398: Output Parameters:
3399: . flg - the flag
3401: Level: developer
3403: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3404: @*/
3405: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3406: {
3408: *flg = s->useFieldOff;
3409: return 0;
3410: }
3412: /*@
3413: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3415: Not Collective
3417: Input Parameters:
3418: + s - the global `PetscSection`
3419: - flg - the flag
3421: Level: developer
3423: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3424: @*/
3425: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3426: {
3428: s->useFieldOff = flg;
3429: return 0;
3430: }
3432: #define PetscSectionExpandPoints_Loop(TYPE) \
3433: { \
3434: PetscInt i, n, o0, o1, size; \
3435: TYPE *a0 = (TYPE *)origArray, *a1; \
3436: PetscSectionGetStorageSize(s, &size); \
3437: PetscMalloc1(size, &a1); \
3438: for (i = 0; i < npoints; i++) { \
3439: PetscSectionGetOffset(origSection, points_[i], &o0); \
3440: PetscSectionGetOffset(s, i, &o1); \
3441: PetscSectionGetDof(s, i, &n); \
3442: PetscMemcpy(&a1[o1], &a0[o0], n *unitsize); \
3443: } \
3444: *newArray = (void *)a1; \
3445: }
3447: /*@
3448: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3450: Not Collective
3452: Input Parameters:
3453: + origSection - the `PetscSection` describing the layout of the array
3454: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
3455: . origArray - the array; its size must be equal to the storage size of origSection
3456: - points - `IS` with points to extract; its indices must lie in the chart of origSection
3458: Output Parameters:
3459: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3460: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3462: Level: developer
3464: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3465: @*/
3466: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3467: {
3468: PetscSection s;
3469: const PetscInt *points_;
3470: PetscInt i, n, npoints, pStart, pEnd;
3471: PetscMPIInt unitsize;
3478: MPI_Type_size(dataType, &unitsize);
3479: ISGetLocalSize(points, &npoints);
3480: ISGetIndices(points, &points_);
3481: PetscSectionGetChart(origSection, &pStart, &pEnd);
3482: PetscSectionCreate(PETSC_COMM_SELF, &s);
3483: PetscSectionSetChart(s, 0, npoints);
3484: for (i = 0; i < npoints; i++) {
3486: PetscSectionGetDof(origSection, points_[i], &n);
3487: PetscSectionSetDof(s, i, n);
3488: }
3489: PetscSectionSetUp(s);
3490: if (newArray) {
3491: if (dataType == MPIU_INT) {
3492: PetscSectionExpandPoints_Loop(PetscInt);
3493: } else if (dataType == MPIU_SCALAR) {
3494: PetscSectionExpandPoints_Loop(PetscScalar);
3495: } else if (dataType == MPIU_REAL) {
3496: PetscSectionExpandPoints_Loop(PetscReal);
3497: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3498: }
3499: if (newSection) {
3500: *newSection = s;
3501: } else {
3502: PetscSectionDestroy(&s);
3503: }
3504: ISRestoreIndices(points, &points_);
3505: return 0;
3506: }