Actual source code: vsectionis.c
petsc-3.9.4 2018-09-11
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/isimpl.h>
6: #include <petscsf.h>
7: #include <petscviewer.h>
9: PetscClassId PETSC_SECTION_CLASSID;
11: /*@
12: PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
14: Collective on MPI_Comm
16: Input Parameters:
17: + comm - the MPI communicator
18: - s - pointer to the section
20: Level: developer
22: Notes:
23: Typical calling sequence
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);
32: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
33: recommended they not be used in user codes unless you really gain something in their use.
35: .seealso: PetscSection, PetscSectionDestroy()
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
43: ISInitializePackage();
45: PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
47: (*s)->pStart = -1;
48: (*s)->pEnd = -1;
49: (*s)->perm = NULL;
50: (*s)->maxDof = 0;
51: (*s)->atlasDof = NULL;
52: (*s)->atlasOff = NULL;
53: (*s)->bc = NULL;
54: (*s)->bcIndices = NULL;
55: (*s)->setup = PETSC_FALSE;
56: (*s)->numFields = 0;
57: (*s)->fieldNames = NULL;
58: (*s)->field = NULL;
59: (*s)->clObj = NULL;
60: (*s)->clSection = NULL;
61: (*s)->clPoints = NULL;
62: (*s)->clSize = 0;
63: (*s)->clPerm = NULL;
64: (*s)->clInvPerm = NULL;
65: return(0);
66: }
68: /*@
69: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
71: Collective on MPI_Comm
73: Input Parameter:
74: . section - the PetscSection
76: Output Parameter:
77: . newSection - the copy
79: Level: developer
81: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
82: @*/
83: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
84: {
85: IS perm;
86: PetscInt numFields, f, pStart, pEnd, p;
92: PetscSectionGetNumFields(section, &numFields);
93: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
94: for (f = 0; f < numFields; ++f) {
95: const char *name = NULL;
96: PetscInt numComp = 0;
98: PetscSectionGetFieldName(section, f, &name);
99: PetscSectionSetFieldName(newSection, f, name);
100: PetscSectionGetFieldComponents(section, f, &numComp);
101: PetscSectionSetFieldComponents(newSection, f, numComp);
102: }
103: PetscSectionGetChart(section, &pStart, &pEnd);
104: PetscSectionSetChart(newSection, pStart, pEnd);
105: PetscSectionGetPermutation(section, &perm);
106: PetscSectionSetPermutation(newSection, perm);
107: for (p = pStart; p < pEnd; ++p) {
108: PetscInt dof, cdof, fcdof = 0;
110: PetscSectionGetDof(section, p, &dof);
111: PetscSectionSetDof(newSection, p, dof);
112: PetscSectionGetConstraintDof(section, p, &cdof);
113: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
114: for (f = 0; f < numFields; ++f) {
115: PetscSectionGetFieldDof(section, p, f, &dof);
116: PetscSectionSetFieldDof(newSection, p, f, dof);
117: if (cdof) {
118: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
119: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
120: }
121: }
122: }
123: PetscSectionSetUp(newSection);
124: for (p = pStart; p < pEnd; ++p) {
125: PetscInt off, cdof, fcdof = 0;
126: const PetscInt *cInd;
128: /* Must set offsets in case they do not agree with the prefix sums */
129: PetscSectionGetOffset(section, p, &off);
130: PetscSectionSetOffset(newSection, p, off);
131: PetscSectionGetConstraintDof(section, p, &cdof);
132: if (cdof) {
133: PetscSectionGetConstraintIndices(section, p, &cInd);
134: PetscSectionSetConstraintIndices(newSection, p, cInd);
135: for (f = 0; f < numFields; ++f) {
136: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
137: if (fcdof) {
138: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
139: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
140: }
141: }
142: }
143: }
144: return(0);
145: }
147: /*@
148: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
150: Collective on MPI_Comm
152: Input Parameter:
153: . section - the PetscSection
155: Output Parameter:
156: . newSection - the copy
158: Level: developer
160: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
161: @*/
162: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
163: {
169: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
170: PetscSectionCopy(section, *newSection);
171: return(0);
172: }
174: /*@
175: PetscSectionCompare - Compares two sections
177: Collective on PetscSection
179: Input Parameterss:
180: + s1 - the first PetscSection
181: - s2 - the second PetscSection
183: Output Parameter:
184: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
186: Level: developer
188: Notes:
189: Field names are disregarded.
191: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
192: @*/
193: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
194: {
195: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
196: const PetscInt *idx1, *idx2;
197: IS perm1, perm2;
198: PetscBool flg;
199: PetscMPIInt mflg;
206: flg = PETSC_FALSE;
208: MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
209: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
210: *congruent = PETSC_FALSE;
211: return(0);
212: }
214: PetscSectionGetChart(s1, &pStart, &pEnd);
215: PetscSectionGetChart(s2, &n1, &n2);
216: if (pStart != n1 || pEnd != n2) goto not_congruent;
218: PetscSectionGetPermutation(s1, &perm1);
219: PetscSectionGetPermutation(s2, &perm2);
220: if (perm1 && perm2) {
221: ISEqual(perm1, perm2, congruent);
222: if (!(*congruent)) goto not_congruent;
223: } else if (perm1 != perm2) goto not_congruent;
225: for (p = pStart; p < pEnd; ++p) {
226: PetscSectionGetOffset(s1, p, &n1);
227: PetscSectionGetOffset(s2, p, &n2);
228: if (n1 != n2) goto not_congruent;
230: PetscSectionGetDof(s1, p, &n1);
231: PetscSectionGetDof(s2, p, &n2);
232: if (n1 != n2) goto not_congruent;
234: PetscSectionGetConstraintDof(s1, p, &ncdof);
235: PetscSectionGetConstraintDof(s2, p, &n2);
236: if (ncdof != n2) goto not_congruent;
238: PetscSectionGetConstraintIndices(s1, p, &idx1);
239: PetscSectionGetConstraintIndices(s2, p, &idx2);
240: PetscMemcmp(idx1, idx2, ncdof*sizeof(PetscInt), congruent);
241: if (!(*congruent)) goto not_congruent;
242: }
244: PetscSectionGetNumFields(s1, &nfields);
245: PetscSectionGetNumFields(s2, &n2);
246: if (nfields != n2) goto not_congruent;
248: for (f = 0; f < nfields; ++f) {
249: PetscSectionGetFieldComponents(s1, f, &n1);
250: PetscSectionGetFieldComponents(s2, f, &n2);
251: if (n1 != n2) goto not_congruent;
253: for (p = pStart; p < pEnd; ++p) {
254: PetscSectionGetFieldOffset(s1, p, f, &n1);
255: PetscSectionGetFieldOffset(s2, p, f, &n2);
256: if (n1 != n2) goto not_congruent;
258: PetscSectionGetFieldDof(s1, p, f, &n1);
259: PetscSectionGetFieldDof(s2, p, f, &n2);
260: if (n1 != n2) goto not_congruent;
262: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
263: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
264: if (nfcdof != n2) goto not_congruent;
266: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
267: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
268: PetscMemcmp(idx1, idx2, nfcdof*sizeof(PetscInt), congruent);
269: if (!(*congruent)) goto not_congruent;
270: }
271: }
273: flg = PETSC_TRUE;
274: not_congruent:
275: MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
276: return(0);
277: }
279: /*@
280: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
282: Not collective
284: Input Parameter:
285: . s - the PetscSection
287: Output Parameter:
288: . numFields - the number of fields defined, or 0 if none were defined
290: Level: intermediate
292: .seealso: PetscSectionSetNumFields()
293: @*/
294: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
295: {
299: *numFields = s->numFields;
300: return(0);
301: }
303: /*@
304: PetscSectionSetNumFields - Sets the number of fields.
306: Not collective
308: Input Parameters:
309: + s - the PetscSection
310: - numFields - the number of fields
312: Level: intermediate
314: .seealso: PetscSectionGetNumFields()
315: @*/
316: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
317: {
318: PetscInt f;
323: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %d must be positive", numFields);
324: PetscSectionReset(s);
326: s->numFields = numFields;
327: PetscMalloc1(s->numFields, &s->numFieldComponents);
328: PetscMalloc1(s->numFields, &s->fieldNames);
329: PetscMalloc1(s->numFields, &s->field);
330: for (f = 0; f < s->numFields; ++f) {
331: char name[64];
333: s->numFieldComponents[f] = 1;
335: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
336: PetscSNPrintf(name, 64, "Field_%D", f);
337: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
338: }
339: return(0);
340: }
342: /*@C
343: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
345: Not Collective
347: Input Parameters:
348: + s - the PetscSection
349: - field - the field number
351: Output Parameter:
352: . fieldName - the field name
354: Level: developer
356: .seealso: PetscSectionSetFieldName()
357: @*/
358: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
359: {
363: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
364: *fieldName = s->fieldNames[field];
365: return(0);
366: }
368: /*@C
369: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
371: Not Collective
373: Input Parameters:
374: + s - the PetscSection
375: . field - the field number
376: - fieldName - the field name
378: Level: developer
380: .seealso: PetscSectionGetFieldName()
381: @*/
382: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
383: {
389: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
390: PetscFree(s->fieldNames[field]);
391: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
392: return(0);
393: }
395: /*@
396: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
398: Not collective
400: Input Parameters:
401: + s - the PetscSection
402: - field - the field number
404: Output Parameter:
405: . numComp - the number of field components
407: Level: intermediate
409: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
410: @*/
411: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
412: {
416: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
417: *numComp = s->numFieldComponents[field];
418: return(0);
419: }
421: /*@
422: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
424: Not collective
426: Input Parameters:
427: + s - the PetscSection
428: . field - the field number
429: - numComp - the number of field components
431: Level: intermediate
433: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
434: @*/
435: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
436: {
439: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
440: s->numFieldComponents[field] = numComp;
441: return(0);
442: }
444: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
445: {
449: if (!s->bc) {
450: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
451: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
452: }
453: return(0);
454: }
456: /*@
457: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
459: Not collective
461: Input Parameter:
462: . s - the PetscSection
464: Output Parameters:
465: + pStart - the first point
466: - pEnd - one past the last point
468: Level: intermediate
470: .seealso: PetscSectionSetChart(), PetscSectionCreate()
471: @*/
472: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
473: {
476: if (pStart) *pStart = s->pStart;
477: if (pEnd) *pEnd = s->pEnd;
478: return(0);
479: }
481: /*@
482: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
484: Not collective
486: Input Parameters:
487: + s - the PetscSection
488: . pStart - the first point
489: - pEnd - one past the last point
491: Level: intermediate
493: .seealso: PetscSectionGetChart(), PetscSectionCreate()
494: @*/
495: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
496: {
497: PetscInt f;
502: /* Cannot Reset() because it destroys field information */
503: s->setup = PETSC_FALSE;
504: PetscSectionDestroy(&s->bc);
505: PetscFree(s->bcIndices);
506: PetscFree2(s->atlasDof, s->atlasOff);
508: s->pStart = pStart;
509: s->pEnd = pEnd;
510: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
511: PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
512: for (f = 0; f < s->numFields; ++f) {
513: PetscSectionSetChart(s->field[f], pStart, pEnd);
514: }
515: return(0);
516: }
518: /*@
519: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
521: Not collective
523: Input Parameter:
524: . s - the PetscSection
526: Output Parameters:
527: . perm - The permutation as an IS
529: Level: intermediate
531: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
532: @*/
533: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
534: {
538: return(0);
539: }
541: /*@
542: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
544: Not collective
546: Input Parameters:
547: + s - the PetscSection
548: - perm - the permutation of points
550: Level: intermediate
552: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
553: @*/
554: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
555: {
561: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
562: if (s->perm != perm) {
563: ISDestroy(&s->perm);
564: if (perm) {
565: s->perm = perm;
566: PetscObjectReference((PetscObject) s->perm);
567: }
568: }
569: return(0);
570: }
572: /*@
573: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
575: Not collective
577: Input Parameters:
578: + s - the PetscSection
579: - point - the point
581: Output Parameter:
582: . numDof - the number of dof
584: Level: intermediate
586: .seealso: PetscSectionSetDof(), PetscSectionCreate()
587: @*/
588: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
589: {
593: #if defined(PETSC_USE_DEBUG)
594: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
595: #endif
596: *numDof = s->atlasDof[point - s->pStart];
597: return(0);
598: }
600: /*@
601: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
603: Not collective
605: Input Parameters:
606: + s - the PetscSection
607: . point - the point
608: - numDof - the number of dof
610: Level: intermediate
612: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
613: @*/
614: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
615: {
618: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
619: s->atlasDof[point - s->pStart] = numDof;
620: return(0);
621: }
623: /*@
624: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
626: Not collective
628: Input Parameters:
629: + s - the PetscSection
630: . point - the point
631: - numDof - the number of additional dof
633: Level: intermediate
635: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
636: @*/
637: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
638: {
641: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
642: s->atlasDof[point - s->pStart] += numDof;
643: return(0);
644: }
646: /*@
647: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
649: Not collective
651: Input Parameters:
652: + s - the PetscSection
653: . point - the point
654: - field - the field
656: Output Parameter:
657: . numDof - the number of dof
659: Level: intermediate
661: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
662: @*/
663: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
664: {
669: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
670: PetscSectionGetDof(s->field[field], point, numDof);
671: return(0);
672: }
674: /*@
675: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
677: Not collective
679: Input Parameters:
680: + s - the PetscSection
681: . point - the point
682: . field - the field
683: - numDof - the number of dof
685: Level: intermediate
687: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
688: @*/
689: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
690: {
695: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
696: PetscSectionSetDof(s->field[field], point, numDof);
697: return(0);
698: }
700: /*@
701: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
703: Not collective
705: Input Parameters:
706: + s - the PetscSection
707: . point - the point
708: . field - the field
709: - numDof - the number of dof
711: Level: intermediate
713: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
714: @*/
715: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
716: {
721: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
722: PetscSectionAddDof(s->field[field], point, numDof);
723: return(0);
724: }
726: /*@
727: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
729: Not collective
731: Input Parameters:
732: + s - the PetscSection
733: - point - the point
735: Output Parameter:
736: . numDof - the number of dof which are fixed by constraints
738: Level: intermediate
740: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
741: @*/
742: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
743: {
749: if (s->bc) {
750: PetscSectionGetDof(s->bc, point, numDof);
751: } else *numDof = 0;
752: return(0);
753: }
755: /*@
756: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
758: Not collective
760: Input Parameters:
761: + s - the PetscSection
762: . point - the point
763: - numDof - the number of dof which are fixed by constraints
765: Level: intermediate
767: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
768: @*/
769: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
770: {
775: if (numDof) {
776: PetscSectionCheckConstraints_Static(s);
777: PetscSectionSetDof(s->bc, point, numDof);
778: }
779: return(0);
780: }
782: /*@
783: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
785: Not collective
787: Input Parameters:
788: + s - the PetscSection
789: . point - the point
790: - numDof - the number of additional dof which are fixed by constraints
792: Level: intermediate
794: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
795: @*/
796: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
797: {
802: if (numDof) {
803: PetscSectionCheckConstraints_Static(s);
804: PetscSectionAddDof(s->bc, point, numDof);
805: }
806: return(0);
807: }
809: /*@
810: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
812: Not collective
814: Input Parameters:
815: + s - the PetscSection
816: . point - the point
817: - field - the field
819: Output Parameter:
820: . numDof - the number of dof which are fixed by constraints
822: Level: intermediate
824: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
825: @*/
826: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
827: {
833: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
834: PetscSectionGetConstraintDof(s->field[field], point, numDof);
835: return(0);
836: }
838: /*@
839: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
841: Not collective
843: Input Parameters:
844: + s - the PetscSection
845: . point - the point
846: . field - the field
847: - numDof - the number of dof which are fixed by constraints
849: Level: intermediate
851: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
852: @*/
853: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
854: {
859: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
860: PetscSectionSetConstraintDof(s->field[field], point, numDof);
861: return(0);
862: }
864: /*@
865: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
867: Not collective
869: Input Parameters:
870: + s - the PetscSection
871: . point - the point
872: . field - the field
873: - numDof - the number of additional dof which are fixed by constraints
875: Level: intermediate
877: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
878: @*/
879: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
880: {
885: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
886: PetscSectionAddConstraintDof(s->field[field], point, numDof);
887: return(0);
888: }
890: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
891: {
896: if (s->bc) {
897: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
899: PetscSectionSetUp(s->bc);
900: PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
901: }
902: return(0);
903: }
905: /*@
906: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
908: Not collective
910: Input Parameter:
911: . s - the PetscSection
913: Level: intermediate
915: .seealso: PetscSectionCreate()
916: @*/
917: PetscErrorCode PetscSectionSetUp(PetscSection s)
918: {
919: const PetscInt *pind = NULL;
920: PetscInt offset = 0, p, f;
921: PetscErrorCode ierr;
925: if (s->setup) return(0);
926: s->setup = PETSC_TRUE;
927: if (s->perm) {ISGetIndices(s->perm, &pind);}
928: for (p = 0; p < s->pEnd - s->pStart; ++p) {
929: const PetscInt q = pind ? pind[p] : p;
931: s->atlasOff[q] = offset;
932: offset += s->atlasDof[q];
933: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
934: }
935: PetscSectionSetUpBC(s);
936: /* Assume that all fields have the same chart */
937: for (p = 0; p < s->pEnd - s->pStart; ++p) {
938: const PetscInt q = pind ? pind[p] : p;
939: PetscInt off = s->atlasOff[q];
941: for (f = 0; f < s->numFields; ++f) {
942: PetscSection sf = s->field[f];
944: sf->atlasOff[q] = off;
945: off += sf->atlasDof[q];
946: }
947: }
948: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
949: for (f = 0; f < s->numFields; ++f) {
950: PetscSectionSetUpBC(s->field[f]);
951: }
952: return(0);
953: }
955: /*@
956: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
958: Not collective
960: Input Parameters:
961: . s - the PetscSection
963: Output Parameter:
964: . maxDof - the maximum dof
966: Level: intermediate
968: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
969: @*/
970: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
971: {
975: *maxDof = s->maxDof;
976: return(0);
977: }
979: /*@
980: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
982: Not collective
984: Input Parameters:
985: + s - the PetscSection
986: - size - the allocated size
988: Output Parameter:
989: . size - the size of an array which can hold all the dofs
991: Level: intermediate
993: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
994: @*/
995: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
996: {
997: PetscInt p, n = 0;
1002: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1003: *size = n;
1004: return(0);
1005: }
1007: /*@
1008: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1010: Not collective
1012: Input Parameters:
1013: + s - the PetscSection
1014: - point - the point
1016: Output Parameter:
1017: . size - the size of an array which can hold all unconstrained dofs
1019: Level: intermediate
1021: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1022: @*/
1023: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1024: {
1025: PetscInt p, n = 0;
1030: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1031: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1032: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1033: }
1034: *size = n;
1035: return(0);
1036: }
1038: /*@
1039: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1040: the local section and an SF describing the section point overlap.
1042: Input Parameters:
1043: + s - The PetscSection for the local field layout
1044: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1045: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1046: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1048: Output Parameter:
1049: . gsection - The PetscSection for the global field layout
1051: Note: This gives negative sizes and offsets to points not owned by this process
1053: Level: developer
1055: .seealso: PetscSectionCreate()
1056: @*/
1057: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1058: {
1059: PetscSection gs;
1060: const PetscInt *pind = NULL;
1061: PetscInt *recv = NULL, *neg = NULL;
1062: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
1063: PetscErrorCode ierr;
1071: PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1072: PetscSectionGetChart(s, &pStart, &pEnd);
1073: PetscSectionSetChart(gs, pStart, pEnd);
1074: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1075: nlocal = nroots; /* The local/leaf space matches global/root space */
1076: /* Must allocate for all points visible to SF, which may be more than this section */
1077: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1078: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
1079: PetscMalloc2(nroots,&neg,nlocal,&recv);
1080: PetscMemzero(neg,nroots*sizeof(PetscInt));
1081: }
1082: /* Mark all local points with negative dof */
1083: for (p = pStart; p < pEnd; ++p) {
1084: PetscSectionGetDof(s, p, &dof);
1085: PetscSectionSetDof(gs, p, dof);
1086: PetscSectionGetConstraintDof(s, p, &cdof);
1087: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1088: if (neg) neg[p] = -(dof+1);
1089: }
1090: PetscSectionSetUpBC(gs);
1091: if (gs->bcIndices) {PetscMemcpy(gs->bcIndices, s->bcIndices, sizeof(PetscInt) * (gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]));}
1092: if (nroots >= 0) {
1093: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1094: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1095: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1096: for (p = pStart; p < pEnd; ++p) {
1097: if (recv[p] < 0) {
1098: gs->atlasDof[p-pStart] = recv[p];
1099: PetscSectionGetDof(s, p, &dof);
1100: if (-(recv[p]+1) != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %d for point %d is not the unconstrained %d", -(recv[p]+1), p, dof);
1101: }
1102: }
1103: }
1104: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1105: if (s->perm) {ISGetIndices(s->perm, &pind);}
1106: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1107: const PetscInt q = pind ? pind[p] : p;
1109: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1110: gs->atlasOff[q] = off;
1111: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1112: }
1113: if (!localOffsets) {
1114: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1115: globalOff -= off;
1116: }
1117: for (p = pStart, off = 0; p < pEnd; ++p) {
1118: gs->atlasOff[p-pStart] += globalOff;
1119: if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1120: }
1121: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1122: /* Put in negative offsets for ghost points */
1123: if (nroots >= 0) {
1124: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1125: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1126: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1127: for (p = pStart; p < pEnd; ++p) {
1128: if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1129: }
1130: }
1131: PetscFree2(neg,recv);
1132: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1133: *gsection = gs;
1134: return(0);
1135: }
1137: /*@
1138: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1139: the local section and an SF describing the section point overlap.
1141: Input Parameters:
1142: + s - The PetscSection for the local field layout
1143: . sf - The SF describing parallel layout of the section points
1144: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1145: . numExcludes - The number of exclusion ranges
1146: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1148: Output Parameter:
1149: . gsection - The PetscSection for the global field layout
1151: Note: This gives negative sizes and offsets to points not owned by this process
1153: Level: developer
1155: .seealso: PetscSectionCreate()
1156: @*/
1157: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1158: {
1159: const PetscInt *pind = NULL;
1160: PetscInt *neg = NULL, *tmpOff = NULL;
1161: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1162: PetscErrorCode ierr;
1168: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1169: PetscSectionGetChart(s, &pStart, &pEnd);
1170: PetscSectionSetChart(*gsection, pStart, pEnd);
1171: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1172: if (nroots >= 0) {
1173: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1174: PetscCalloc1(nroots, &neg);
1175: if (nroots > pEnd-pStart) {
1176: PetscCalloc1(nroots, &tmpOff);
1177: } else {
1178: tmpOff = &(*gsection)->atlasDof[-pStart];
1179: }
1180: }
1181: /* Mark ghost points with negative dof */
1182: for (p = pStart; p < pEnd; ++p) {
1183: for (e = 0; e < numExcludes; ++e) {
1184: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1185: PetscSectionSetDof(*gsection, p, 0);
1186: break;
1187: }
1188: }
1189: if (e < numExcludes) continue;
1190: PetscSectionGetDof(s, p, &dof);
1191: PetscSectionSetDof(*gsection, p, dof);
1192: PetscSectionGetConstraintDof(s, p, &cdof);
1193: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1194: if (neg) neg[p] = -(dof+1);
1195: }
1196: PetscSectionSetUpBC(*gsection);
1197: if (nroots >= 0) {
1198: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1199: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1200: if (nroots > pEnd - pStart) {
1201: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1202: }
1203: }
1204: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1205: if (s->perm) {ISGetIndices(s->perm, &pind);}
1206: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1207: const PetscInt q = pind ? pind[p] : p;
1209: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1210: (*gsection)->atlasOff[q] = off;
1211: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1212: }
1213: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1214: globalOff -= off;
1215: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1216: (*gsection)->atlasOff[p] += globalOff;
1217: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1218: }
1219: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1220: /* Put in negative offsets for ghost points */
1221: if (nroots >= 0) {
1222: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1223: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1224: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1225: if (nroots > pEnd - pStart) {
1226: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1227: }
1228: }
1229: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1230: PetscFree(neg);
1231: return(0);
1232: }
1234: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1235: {
1236: PetscInt pStart, pEnd, p, localSize = 0;
1240: PetscSectionGetChart(s, &pStart, &pEnd);
1241: for (p = pStart; p < pEnd; ++p) {
1242: PetscInt dof;
1244: PetscSectionGetDof(s, p, &dof);
1245: if (dof > 0) ++localSize;
1246: }
1247: PetscLayoutCreate(comm, layout);
1248: PetscLayoutSetLocalSize(*layout, localSize);
1249: PetscLayoutSetBlockSize(*layout, 1);
1250: PetscLayoutSetUp(*layout);
1251: return(0);
1252: }
1254: /*@
1255: PetscSectionGetValueLayout - Get the PetscLayout associated with a section, usually the default global section.
1257: Input Parameters:
1258: + comm - The MPI_Comm
1259: - s - The PetscSection
1261: Output Parameter:
1262: . layout - The layout for the section
1264: Level: developer
1266: .seealso: PetscSectionCreate()
1267: @*/
1268: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1269: {
1270: PetscInt pStart, pEnd, p, localSize = 0;
1276: PetscSectionGetChart(s, &pStart, &pEnd);
1277: for (p = pStart; p < pEnd; ++p) {
1278: PetscInt dof,cdof;
1280: PetscSectionGetDof(s, p, &dof);
1281: PetscSectionGetConstraintDof(s, p, &cdof);
1282: if (dof-cdof > 0) localSize += dof-cdof;
1283: }
1284: PetscLayoutCreate(comm, layout);
1285: PetscLayoutSetLocalSize(*layout, localSize);
1286: PetscLayoutSetBlockSize(*layout, 1);
1287: PetscLayoutSetUp(*layout);
1288: return(0);
1289: }
1291: /*@
1292: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1294: Not collective
1296: Input Parameters:
1297: + s - the PetscSection
1298: - point - the point
1300: Output Parameter:
1301: . offset - the offset
1303: Level: intermediate
1305: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1306: @*/
1307: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1308: {
1312: #if defined(PETSC_USE_DEBUG)
1313: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
1314: #endif
1315: *offset = s->atlasOff[point - s->pStart];
1316: return(0);
1317: }
1319: /*@
1320: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1322: Not collective
1324: Input Parameters:
1325: + s - the PetscSection
1326: . point - the point
1327: - offset - the offset
1329: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1331: Level: intermediate
1333: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1334: @*/
1335: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1336: {
1339: if ((point < s->pStart) || (point >= s->pEnd)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %d should be in [%d, %d)", point, s->pStart, s->pEnd);
1340: s->atlasOff[point - s->pStart] = offset;
1341: return(0);
1342: }
1344: /*@
1345: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1347: Not collective
1349: Input Parameters:
1350: + s - the PetscSection
1351: . point - the point
1352: - field - the field
1354: Output Parameter:
1355: . offset - the offset
1357: Level: intermediate
1359: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1360: @*/
1361: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1362: {
1368: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1369: PetscSectionGetOffset(s->field[field], point, offset);
1370: return(0);
1371: }
1373: /*@
1374: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1376: Not collective
1378: Input Parameters:
1379: + s - the PetscSection
1380: . point - the point
1381: . field - the field
1382: - offset - the offset
1384: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1386: Level: intermediate
1388: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1389: @*/
1390: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1391: {
1396: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1397: PetscSectionSetOffset(s->field[field], point, offset);
1398: return(0);
1399: }
1401: /* This gives the offset on a point of the field, ignoring constraints */
1402: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1403: {
1404: PetscInt off, foff;
1410: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
1411: PetscSectionGetOffset(s, point, &off);
1412: PetscSectionGetOffset(s->field[field], point, &foff);
1413: *offset = foff - off;
1414: return(0);
1415: }
1417: /*@
1418: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1420: Not collective
1422: Input Parameter:
1423: . s - the PetscSection
1425: Output Parameters:
1426: + start - the minimum offset
1427: - end - one more than the maximum offset
1429: Level: intermediate
1431: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1432: @*/
1433: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1434: {
1435: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1440: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1441: PetscSectionGetChart(s, &pStart, &pEnd);
1442: for (p = 0; p < pEnd-pStart; ++p) {
1443: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1445: if (off >= 0) {
1446: os = PetscMin(os, off);
1447: oe = PetscMax(oe, off+dof);
1448: }
1449: }
1450: if (start) *start = os;
1451: if (end) *end = oe;
1452: return(0);
1453: }
1455: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, const PetscInt fields[], PetscSection *subs)
1456: {
1457: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1461: if (!numFields) return(0);
1465: PetscSectionGetNumFields(s, &nF);
1466: if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1467: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1468: PetscSectionSetNumFields(*subs, numFields);
1469: for (f = 0; f < numFields; ++f) {
1470: const char *name = NULL;
1471: PetscInt numComp = 0;
1473: PetscSectionGetFieldName(s, fields[f], &name);
1474: PetscSectionSetFieldName(*subs, f, name);
1475: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1476: PetscSectionSetFieldComponents(*subs, f, numComp);
1477: }
1478: PetscSectionGetChart(s, &pStart, &pEnd);
1479: PetscSectionSetChart(*subs, pStart, pEnd);
1480: for (p = pStart; p < pEnd; ++p) {
1481: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1483: for (f = 0; f < numFields; ++f) {
1484: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1485: PetscSectionSetFieldDof(*subs, p, f, fdof);
1486: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1487: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1488: dof += fdof;
1489: cdof += cfdof;
1490: }
1491: PetscSectionSetDof(*subs, p, dof);
1492: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1493: maxCdof = PetscMax(cdof, maxCdof);
1494: }
1495: PetscSectionSetUp(*subs);
1496: if (maxCdof) {
1497: PetscInt *indices;
1499: PetscMalloc1(maxCdof, &indices);
1500: for (p = pStart; p < pEnd; ++p) {
1501: PetscInt cdof;
1503: PetscSectionGetConstraintDof(*subs, p, &cdof);
1504: if (cdof) {
1505: const PetscInt *oldIndices = NULL;
1506: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1508: for (f = 0; f < numFields; ++f) {
1509: PetscInt oldFoff = 0, oldf;
1511: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1512: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1513: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1514: /* This can be sped up if we assume sorted fields */
1515: for (oldf = 0; oldf < fields[f]; ++oldf) {
1516: PetscInt oldfdof = 0;
1517: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1518: oldFoff += oldfdof;
1519: }
1520: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1521: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1522: numConst += cfdof;
1523: fOff += fdof;
1524: }
1525: PetscSectionSetConstraintIndices(*subs, p, indices);
1526: }
1527: }
1528: PetscFree(indices);
1529: }
1530: return(0);
1531: }
1533: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1534: {
1535: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1539: if (!len) return(0);
1540: for (i = 0; i < len; ++i) {
1541: PetscInt nf, pStarti, pEndi;
1543: PetscSectionGetNumFields(s[i], &nf);
1544: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1545: pStart = PetscMin(pStart, pStarti);
1546: pEnd = PetscMax(pEnd, pEndi);
1547: Nf += nf;
1548: }
1549: PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1550: PetscSectionSetNumFields(*supers, Nf);
1551: for (i = 0, f = 0; i < len; ++i) {
1552: PetscInt nf, fi;
1554: PetscSectionGetNumFields(s[i], &nf);
1555: for (fi = 0; fi < nf; ++fi, ++f) {
1556: const char *name = NULL;
1557: PetscInt numComp = 0;
1559: PetscSectionGetFieldName(s[i], fi, &name);
1560: PetscSectionSetFieldName(*supers, f, name);
1561: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1562: PetscSectionSetFieldComponents(*supers, f, numComp);
1563: }
1564: }
1565: PetscSectionSetChart(*supers, pStart, pEnd);
1566: for (p = pStart; p < pEnd; ++p) {
1567: PetscInt dof = 0, cdof = 0;
1569: for (i = 0, f = 0; i < len; ++i) {
1570: PetscInt nf, fi, pStarti, pEndi;
1571: PetscInt fdof = 0, cfdof = 0;
1573: PetscSectionGetNumFields(s[i], &nf);
1574: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1575: if ((p < pStarti) || (p >= pEndi)) continue;
1576: for (fi = 0; fi < nf; ++fi, ++f) {
1577: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1578: PetscSectionAddFieldDof(*supers, p, f, fdof);
1579: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1580: if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1581: dof += fdof;
1582: cdof += cfdof;
1583: }
1584: }
1585: PetscSectionSetDof(*supers, p, dof);
1586: if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1587: maxCdof = PetscMax(cdof, maxCdof);
1588: }
1589: PetscSectionSetUp(*supers);
1590: if (maxCdof) {
1591: PetscInt *indices;
1593: PetscMalloc1(maxCdof, &indices);
1594: for (p = pStart; p < pEnd; ++p) {
1595: PetscInt cdof;
1597: PetscSectionGetConstraintDof(*supers, p, &cdof);
1598: if (cdof) {
1599: PetscInt dof, numConst = 0, fOff = 0;
1601: for (i = 0, f = 0; i < len; ++i) {
1602: const PetscInt *oldIndices = NULL;
1603: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1605: PetscSectionGetNumFields(s[i], &nf);
1606: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1607: if ((p < pStarti) || (p >= pEndi)) continue;
1608: for (fi = 0; fi < nf; ++fi, ++f) {
1609: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1610: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1611: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1612: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1613: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1614: numConst += cfdof;
1615: }
1616: PetscSectionGetDof(s[i], p, &dof);
1617: fOff += dof;
1618: }
1619: PetscSectionSetConstraintIndices(*supers, p, indices);
1620: }
1621: }
1622: PetscFree(indices);
1623: }
1624: return(0);
1625: }
1627: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1628: {
1629: const PetscInt *points = NULL, *indices = NULL;
1630: PetscInt numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;
1637: PetscSectionGetNumFields(s, &numFields);
1638: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1639: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1640: for (f = 0; f < numFields; ++f) {
1641: const char *name = NULL;
1642: PetscInt numComp = 0;
1644: PetscSectionGetFieldName(s, f, &name);
1645: PetscSectionSetFieldName(*subs, f, name);
1646: PetscSectionGetFieldComponents(s, f, &numComp);
1647: PetscSectionSetFieldComponents(*subs, f, numComp);
1648: }
1649: /* For right now, we do not try to squeeze the subchart */
1650: if (subpointMap) {
1651: ISGetSize(subpointMap, &numSubpoints);
1652: ISGetIndices(subpointMap, &points);
1653: }
1654: PetscSectionGetChart(s, &pStart, &pEnd);
1655: PetscSectionSetChart(*subs, 0, numSubpoints);
1656: for (p = pStart; p < pEnd; ++p) {
1657: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1659: PetscFindInt(p, numSubpoints, points, &subp);
1660: if (subp < 0) continue;
1661: for (f = 0; f < numFields; ++f) {
1662: PetscSectionGetFieldDof(s, p, f, &fdof);
1663: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1664: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1665: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1666: }
1667: PetscSectionGetDof(s, p, &dof);
1668: PetscSectionSetDof(*subs, subp, dof);
1669: PetscSectionGetConstraintDof(s, p, &cdof);
1670: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1671: }
1672: PetscSectionSetUp(*subs);
1673: /* Change offsets to original offsets */
1674: for (p = pStart; p < pEnd; ++p) {
1675: PetscInt off, foff = 0;
1677: PetscFindInt(p, numSubpoints, points, &subp);
1678: if (subp < 0) continue;
1679: for (f = 0; f < numFields; ++f) {
1680: PetscSectionGetFieldOffset(s, p, f, &foff);
1681: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1682: }
1683: PetscSectionGetOffset(s, p, &off);
1684: PetscSectionSetOffset(*subs, subp, off);
1685: }
1686: /* Copy constraint indices */
1687: for (subp = 0; subp < numSubpoints; ++subp) {
1688: PetscInt cdof;
1690: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1691: if (cdof) {
1692: for (f = 0; f < numFields; ++f) {
1693: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1694: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1695: }
1696: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1697: PetscSectionSetConstraintIndices(*subs, subp, indices);
1698: }
1699: }
1700: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1701: return(0);
1702: }
1704: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1705: {
1706: PetscInt p;
1707: PetscMPIInt rank;
1711: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1712: PetscViewerASCIIPushSynchronized(viewer);
1713: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1714: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1715: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1716: PetscInt b;
1718: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1719: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1720: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1721: }
1722: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1723: } else {
1724: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1725: }
1726: }
1727: PetscViewerFlush(viewer);
1728: PetscViewerASCIIPopSynchronized(viewer);
1729: if (s->sym) {
1730: PetscViewerASCIIPushTab(viewer);
1731: PetscSectionSymView(s->sym,viewer);
1732: PetscViewerASCIIPopTab(viewer);
1733: }
1734: return(0);
1735: }
1737: /*@C
1738: PetscSectionView - Views a PetscSection
1740: Collective on PetscSection
1742: Input Parameters:
1743: + s - the PetscSection object to view
1744: - v - the viewer
1746: Level: developer
1748: .seealso PetscSectionCreate(), PetscSectionDestroy()
1749: @*/
1750: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1751: {
1752: PetscBool isascii;
1753: PetscInt f;
1758: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1760: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1761: if (isascii) {
1762: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1763: if (s->numFields) {
1764: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1765: for (f = 0; f < s->numFields; ++f) {
1766: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
1767: PetscSectionView_ASCII(s->field[f], viewer);
1768: }
1769: } else {
1770: PetscSectionView_ASCII(s, viewer);
1771: }
1772: }
1773: return(0);
1774: }
1776: /*@
1777: PetscSectionReset - Frees all section data.
1779: Not collective
1781: Input Parameters:
1782: . s - the PetscSection
1784: Level: developer
1786: .seealso: PetscSection, PetscSectionCreate()
1787: @*/
1788: PetscErrorCode PetscSectionReset(PetscSection s)
1789: {
1790: PetscInt f;
1795: PetscFree(s->numFieldComponents);
1796: for (f = 0; f < s->numFields; ++f) {
1797: PetscSectionDestroy(&s->field[f]);
1798: PetscFree(s->fieldNames[f]);
1799: }
1800: PetscFree(s->fieldNames);
1801: PetscFree(s->field);
1802: PetscSectionDestroy(&s->bc);
1803: PetscFree(s->bcIndices);
1804: PetscFree2(s->atlasDof, s->atlasOff);
1805: PetscSectionDestroy(&s->clSection);
1806: ISDestroy(&s->clPoints);
1807: ISDestroy(&s->perm);
1808: PetscFree(s->clPerm);
1809: PetscFree(s->clInvPerm);
1810: PetscSectionSymDestroy(&s->sym);
1812: s->pStart = -1;
1813: s->pEnd = -1;
1814: s->maxDof = 0;
1815: s->setup = PETSC_FALSE;
1816: s->numFields = 0;
1817: s->clObj = NULL;
1818: return(0);
1819: }
1821: /*@
1822: PetscSectionDestroy - Frees a section object and frees its range if that exists.
1824: Not collective
1826: Input Parameters:
1827: . s - the PetscSection
1829: Level: developer
1831: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1832: recommended they not be used in user codes unless you really gain something in their use.
1834: .seealso: PetscSection, PetscSectionCreate()
1835: @*/
1836: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1837: {
1841: if (!*s) return(0);
1843: if (--((PetscObject)(*s))->refct > 0) {
1844: *s = NULL;
1845: return(0);
1846: }
1847: PetscSectionReset(*s);
1848: PetscHeaderDestroy(s);
1849: return(0);
1850: }
1852: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1853: {
1854: const PetscInt p = point - s->pStart;
1858: *values = &baseArray[s->atlasOff[p]];
1859: return(0);
1860: }
1862: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1863: {
1864: PetscInt *array;
1865: const PetscInt p = point - s->pStart;
1866: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1867: PetscInt cDim = 0;
1872: PetscSectionGetConstraintDof(s, p, &cDim);
1873: array = &baseArray[s->atlasOff[p]];
1874: if (!cDim) {
1875: if (orientation >= 0) {
1876: const PetscInt dim = s->atlasDof[p];
1877: PetscInt i;
1879: if (mode == INSERT_VALUES) {
1880: for (i = 0; i < dim; ++i) array[i] = values[i];
1881: } else {
1882: for (i = 0; i < dim; ++i) array[i] += values[i];
1883: }
1884: } else {
1885: PetscInt offset = 0;
1886: PetscInt j = -1, field, i;
1888: for (field = 0; field < s->numFields; ++field) {
1889: const PetscInt dim = s->field[field]->atlasDof[p];
1891: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1892: offset += dim;
1893: }
1894: }
1895: } else {
1896: if (orientation >= 0) {
1897: const PetscInt dim = s->atlasDof[p];
1898: PetscInt cInd = 0, i;
1899: const PetscInt *cDof;
1901: PetscSectionGetConstraintIndices(s, point, &cDof);
1902: if (mode == INSERT_VALUES) {
1903: for (i = 0; i < dim; ++i) {
1904: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1905: array[i] = values[i];
1906: }
1907: } else {
1908: for (i = 0; i < dim; ++i) {
1909: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1910: array[i] += values[i];
1911: }
1912: }
1913: } else {
1914: const PetscInt *cDof;
1915: PetscInt offset = 0;
1916: PetscInt cOffset = 0;
1917: PetscInt j = 0, field;
1919: PetscSectionGetConstraintIndices(s, point, &cDof);
1920: for (field = 0; field < s->numFields; ++field) {
1921: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
1922: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1923: const PetscInt sDim = dim - tDim;
1924: PetscInt cInd = 0, i,k;
1926: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1927: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1928: array[j] = values[k];
1929: }
1930: offset += dim;
1931: cOffset += dim - tDim;
1932: }
1933: }
1934: }
1935: return(0);
1936: }
1938: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1939: {
1943: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1944: return(0);
1945: }
1947: /*@C
1948: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
1950: Input Parameters:
1951: + s - The PetscSection
1952: - point - The point
1954: Output Parameter:
1955: . indices - The constrained dofs
1957: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
1959: Level: advanced
1961: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1962: @*/
1963: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1964: {
1969: if (s->bc) {
1970: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1971: } else *indices = NULL;
1972: return(0);
1973: }
1975: /*@C
1976: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
1978: Input Parameters:
1979: + s - The PetscSection
1980: . point - The point
1981: - indices - The constrained dofs
1983: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
1985: Level: advanced
1987: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1988: @*/
1989: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1990: {
1995: if (s->bc) {
1996: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1997: }
1998: return(0);
1999: }
2001: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2002: {
2007: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
2008: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2009: return(0);
2010: }
2012: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2013: {
2018: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
2019: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2020: return(0);
2021: }
2023: /*@
2024: PetscSectionPermute - Reorder the section according to the input point permutation
2026: Collective on PetscSection
2028: Input Parameter:
2029: + section - The PetscSection object
2030: - perm - The point permutation, old point p becomes new point perm[p]
2032: Output Parameter:
2033: . sectionNew - The permuted PetscSection
2035: Level: intermediate
2037: .keywords: mesh
2038: .seealso: MatPermute()
2039: @*/
2040: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2041: {
2042: PetscSection s = section, sNew;
2043: const PetscInt *perm;
2044: PetscInt numFields, f, numPoints, pStart, pEnd, p;
2045: PetscErrorCode ierr;
2051: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2052: PetscSectionGetNumFields(s, &numFields);
2053: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2054: for (f = 0; f < numFields; ++f) {
2055: const char *name;
2056: PetscInt numComp;
2058: PetscSectionGetFieldName(s, f, &name);
2059: PetscSectionSetFieldName(sNew, f, name);
2060: PetscSectionGetFieldComponents(s, f, &numComp);
2061: PetscSectionSetFieldComponents(sNew, f, numComp);
2062: }
2063: ISGetLocalSize(permutation, &numPoints);
2064: ISGetIndices(permutation, &perm);
2065: PetscSectionGetChart(s, &pStart, &pEnd);
2066: PetscSectionSetChart(sNew, pStart, pEnd);
2067: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
2068: for (p = pStart; p < pEnd; ++p) {
2069: PetscInt dof, cdof;
2071: PetscSectionGetDof(s, p, &dof);
2072: PetscSectionSetDof(sNew, perm[p], dof);
2073: PetscSectionGetConstraintDof(s, p, &cdof);
2074: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2075: for (f = 0; f < numFields; ++f) {
2076: PetscSectionGetFieldDof(s, p, f, &dof);
2077: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2078: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2079: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2080: }
2081: }
2082: PetscSectionSetUp(sNew);
2083: for (p = pStart; p < pEnd; ++p) {
2084: const PetscInt *cind;
2085: PetscInt cdof;
2087: PetscSectionGetConstraintDof(s, p, &cdof);
2088: if (cdof) {
2089: PetscSectionGetConstraintIndices(s, p, &cind);
2090: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2091: }
2092: for (f = 0; f < numFields; ++f) {
2093: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2094: if (cdof) {
2095: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2096: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2097: }
2098: }
2099: }
2100: ISRestoreIndices(permutation, &perm);
2101: *sectionNew = sNew;
2102: return(0);
2103: }
2105: /*@C
2106: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
2108: Collective
2110: Input Parameters:
2111: + sf - The SF
2112: - rootSection - Section defined on root space
2114: Output Parameters:
2115: + remoteOffsets - root offsets in leaf storage, or NULL
2116: - leafSection - Section defined on the leaf space
2118: Level: intermediate
2120: .seealso: PetscSFCreate()
2121: @*/
2122: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2123: {
2124: PetscSF embedSF;
2125: const PetscInt *ilocal, *indices;
2126: IS selected;
2127: PetscInt numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;
2131: PetscSectionGetNumFields(rootSection, &numFields);
2132: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2133: for (f = 0; f < numFields; ++f) {
2134: PetscInt numComp = 0;
2135: PetscSectionGetFieldComponents(rootSection, f, &numComp);
2136: PetscSectionSetFieldComponents(leafSection, f, numComp);
2137: }
2138: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2139: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2140: rpEnd = PetscMin(rpEnd,nroots);
2141: rpEnd = PetscMax(rpStart,rpEnd);
2142: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2143: ISGetIndices(selected, &indices);
2144: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2145: ISRestoreIndices(selected, &indices);
2146: ISDestroy(&selected);
2147: PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
2148: if (nleaves && ilocal) {
2149: for (i = 0; i < nleaves; ++i) {
2150: lpStart = PetscMin(lpStart, ilocal[i]);
2151: lpEnd = PetscMax(lpEnd, ilocal[i]);
2152: }
2153: ++lpEnd;
2154: } else {
2155: lpStart = 0;
2156: lpEnd = nleaves;
2157: }
2158: PetscSectionSetChart(leafSection, lpStart, lpEnd);
2159: /* Could fuse these at the cost of a copy and extra allocation */
2160: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2161: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2162: for (f = 0; f < numFields; ++f) {
2163: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2164: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2165: }
2166: if (remoteOffsets) {
2167: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2168: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2169: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2170: }
2171: PetscSFDestroy(&embedSF);
2172: PetscSectionSetUp(leafSection);
2173: return(0);
2174: }
2176: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2177: {
2178: PetscSF embedSF;
2179: const PetscInt *indices;
2180: IS selected;
2181: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2182: PetscErrorCode ierr;
2185: *remoteOffsets = NULL;
2186: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2187: if (numRoots < 0) return(0);
2188: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2189: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2190: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2191: ISGetIndices(selected, &indices);
2192: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2193: ISRestoreIndices(selected, &indices);
2194: ISDestroy(&selected);
2195: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2196: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2197: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2198: PetscSFDestroy(&embedSF);
2199: return(0);
2200: }
2202: /*@C
2203: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2205: Input Parameters:
2206: + sf - The SF
2207: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2208: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2209: - leafSection - Data layout of local points for incoming data (this is the distributed section)
2211: Output Parameters:
2212: - sectionSF - The new SF
2214: Note: Either rootSection or remoteOffsets can be specified
2216: Level: intermediate
2218: .seealso: PetscSFCreate()
2219: @*/
2220: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2221: {
2222: MPI_Comm comm;
2223: const PetscInt *localPoints;
2224: const PetscSFNode *remotePoints;
2225: PetscInt lpStart, lpEnd;
2226: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2227: PetscInt *localIndices;
2228: PetscSFNode *remoteIndices;
2229: PetscInt i, ind;
2230: PetscErrorCode ierr;
2238: PetscObjectGetComm((PetscObject)sf,&comm);
2239: PetscSFCreate(comm, sectionSF);
2240: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2241: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2242: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2243: if (numRoots < 0) return(0);
2244: for (i = 0; i < numPoints; ++i) {
2245: PetscInt localPoint = localPoints ? localPoints[i] : i;
2246: PetscInt dof;
2248: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2249: PetscSectionGetDof(leafSection, localPoint, &dof);
2250: numIndices += dof;
2251: }
2252: }
2253: PetscMalloc1(numIndices, &localIndices);
2254: PetscMalloc1(numIndices, &remoteIndices);
2255: /* Create new index graph */
2256: for (i = 0, ind = 0; i < numPoints; ++i) {
2257: PetscInt localPoint = localPoints ? localPoints[i] : i;
2258: PetscInt rank = remotePoints[i].rank;
2260: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2261: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2262: PetscInt loff, dof, d;
2264: PetscSectionGetOffset(leafSection, localPoint, &loff);
2265: PetscSectionGetDof(leafSection, localPoint, &dof);
2266: for (d = 0; d < dof; ++d, ++ind) {
2267: localIndices[ind] = loff+d;
2268: remoteIndices[ind].rank = rank;
2269: remoteIndices[ind].index = remoteOffset+d;
2270: }
2271: }
2272: }
2273: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2274: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2275: return(0);
2276: }
2278: /*@
2279: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2281: Input Parameters:
2282: + section - The PetscSection
2283: . obj - A PetscObject which serves as the key for this index
2284: . clSection - Section giving the size of the closure of each point
2285: - clPoints - IS giving the points in each closure
2287: Note: We compress out closure points with no dofs in this section
2289: Level: intermediate
2291: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2292: @*/
2293: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2294: {
2298: if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2299: section->clObj = obj;
2300: PetscSectionDestroy(§ion->clSection);
2301: ISDestroy(§ion->clPoints);
2302: section->clSection = clSection;
2303: section->clPoints = clPoints;
2304: return(0);
2305: }
2307: /*@
2308: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2310: Input Parameters:
2311: + section - The PetscSection
2312: - obj - A PetscObject which serves as the key for this index
2314: Output Parameters:
2315: + clSection - Section giving the size of the closure of each point
2316: - clPoints - IS giving the points in each closure
2318: Note: We compress out closure points with no dofs in this section
2320: Level: intermediate
2322: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2323: @*/
2324: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2325: {
2327: if (section->clObj == obj) {
2328: if (clSection) *clSection = section->clSection;
2329: if (clPoints) *clPoints = section->clPoints;
2330: } else {
2331: if (clSection) *clSection = NULL;
2332: if (clPoints) *clPoints = NULL;
2333: }
2334: return(0);
2335: }
2337: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2338: {
2339: PetscInt i;
2343: if (section->clObj != obj) {
2344: PetscSectionDestroy(§ion->clSection);
2345: ISDestroy(§ion->clPoints);
2346: }
2347: section->clObj = obj;
2348: PetscFree(section->clPerm);
2349: PetscFree(section->clInvPerm);
2350: section->clSize = clSize;
2351: if (mode == PETSC_COPY_VALUES) {
2352: PetscMalloc1(clSize, §ion->clPerm);
2353: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2354: PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2355: } else if (mode == PETSC_OWN_POINTER) {
2356: section->clPerm = clPerm;
2357: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2358: PetscMalloc1(clSize, §ion->clInvPerm);
2359: for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2360: return(0);
2361: }
2363: /*@
2364: PetscSectionSetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2366: Input Parameters:
2367: + section - The PetscSection
2368: . obj - A PetscObject which serves as the key for this index
2369: - perm - Permutation of the cell dof closure
2371: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2372: other points (like faces).
2374: Level: intermediate
2376: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2377: @*/
2378: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2379: {
2380: const PetscInt *clPerm = NULL;
2381: PetscInt clSize = 0;
2382: PetscErrorCode ierr;
2385: if (perm) {
2386: ISGetLocalSize(perm, &clSize);
2387: ISGetIndices(perm, &clPerm);
2388: }
2389: PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2390: if (perm) {ISRestoreIndices(perm, &clPerm);}
2391: return(0);
2392: }
2394: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2395: {
2397: if (section->clObj == obj) {
2398: if (size) *size = section->clSize;
2399: if (perm) *perm = section->clPerm;
2400: } else {
2401: if (size) *size = 0;
2402: if (perm) *perm = NULL;
2403: }
2404: return(0);
2405: }
2407: /*@
2408: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2410: Input Parameters:
2411: + section - The PetscSection
2412: - obj - A PetscObject which serves as the key for this index
2414: Output Parameter:
2415: . perm - The dof closure permutation
2417: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2418: other points (like faces).
2420: The user must destroy the IS that is returned.
2422: Level: intermediate
2424: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2425: @*/
2426: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2427: {
2428: const PetscInt *clPerm;
2429: PetscInt clSize;
2430: PetscErrorCode ierr;
2433: PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2434: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2435: return(0);
2436: }
2438: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2439: {
2441: if (section->clObj == obj) {
2442: if (size) *size = section->clSize;
2443: if (perm) *perm = section->clInvPerm;
2444: } else {
2445: if (size) *size = 0;
2446: if (perm) *perm = NULL;
2447: }
2448: return(0);
2449: }
2451: /*@
2452: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2454: Input Parameters:
2455: + section - The PetscSection
2456: - obj - A PetscObject which serves as the key for this index
2458: Output Parameters:
2459: + size - The dof closure size
2460: - perm - The dof closure permutation
2462: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2463: other points (like faces).
2465: The user must destroy the IS that is returned.
2467: Level: intermediate
2469: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2470: @*/
2471: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2472: {
2473: const PetscInt *clPerm;
2474: PetscInt clSize;
2475: PetscErrorCode ierr;
2478: PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2479: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2480: return(0);
2481: }
2483: /*@
2484: PetscSectionGetField - Get the subsection associated with a single field
2486: Input Parameters:
2487: + s - The PetscSection
2488: - field - The field number
2490: Output Parameter:
2491: . subs - The subsection for the given field
2493: Level: intermediate
2495: .seealso: PetscSectionSetNumFields()
2496: @*/
2497: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2498: {
2504: if ((field < 0) || (field >= s->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field %d should be in [%d, %d)", field, 0, s->numFields);
2505: PetscObjectReference((PetscObject) s->field[field]);
2506: *subs = s->field[field];
2507: return(0);
2508: }
2510: PetscClassId PETSC_SECTION_SYM_CLASSID;
2511: PetscFunctionList PetscSectionSymList = NULL;
2513: /*@
2514: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
2516: Collective on MPI_Comm
2518: Input Parameter:
2519: . comm - the MPI communicator
2521: Output Parameter:
2522: . sym - pointer to the new set of symmetries
2524: Level: developer
2526: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2527: @*/
2528: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2529: {
2534: ISInitializePackage();
2535: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2536: return(0);
2537: }
2539: /*@C
2540: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
2542: Collective on PetscSectionSym
2544: Input Parameters:
2545: + sym - The section symmetry object
2546: - method - The name of the section symmetry type
2548: Level: developer
2550: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2551: @*/
2552: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2553: {
2554: PetscErrorCode (*r)(PetscSectionSym);
2555: PetscBool match;
2560: PetscObjectTypeCompare((PetscObject) sym, method, &match);
2561: if (match) return(0);
2563: PetscFunctionListFind(PetscSectionSymList,method,&r);
2564: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2565: if (sym->ops->destroy) {
2566: (*sym->ops->destroy)(sym);
2567: sym->ops->destroy = NULL;
2568: }
2569: (*r)(sym);
2570: PetscObjectChangeTypeName((PetscObject)sym,method);
2571: return(0);
2572: }
2575: /*@C
2576: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
2578: Not Collective
2580: Input Parameter:
2581: . sym - The section symmetry
2583: Output Parameter:
2584: . type - The index set type name
2586: Level: developer
2588: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2589: @*/
2590: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2591: {
2595: *type = ((PetscObject)sym)->type_name;
2596: return(0);
2597: }
2599: /*@C
2600: PetscSectionSymRegister - Adds a new section symmetry implementation
2602: Not Collective
2604: Input Parameters:
2605: + name - The name of a new user-defined creation routine
2606: - create_func - The creation routine itself
2608: Notes:
2609: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
2611: Level: developer
2613: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2614: @*/
2615: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2616: {
2620: ISInitializePackage();
2621: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2622: return(0);
2623: }
2625: /*@
2626: PetscSectionSymDestroy - Destroys a section symmetry.
2628: Collective on PetscSectionSym
2630: Input Parameters:
2631: . sym - the section symmetry
2633: Level: developer
2635: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2636: @*/
2637: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2638: {
2639: SymWorkLink link,next;
2643: if (!*sym) return(0);
2645: if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2646: if ((*sym)->ops->destroy) {
2647: (*(*sym)->ops->destroy)(*sym);
2648: }
2649: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2650: for (link=(*sym)->workin; link; link=next) {
2651: next = link->next;
2652: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2653: PetscFree(link);
2654: }
2655: (*sym)->workin = NULL;
2656: PetscHeaderDestroy(sym);
2657: return(0);
2658: }
2660: /*@C
2661: PetscSectionSymView - Displays a section symmetry
2663: Collective on PetscSectionSym
2665: Input Parameters:
2666: + sym - the index set
2667: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
2669: Level: developer
2671: .seealso: PetscViewerASCIIOpen()
2672: @*/
2673: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2674: {
2679: if (!viewer) {
2680: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2681: }
2684: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2685: if (sym->ops->view) {
2686: (*sym->ops->view)(sym,viewer);
2687: }
2688: return(0);
2689: }
2691: /*@
2692: PetscSectionSetSym - Set the symmetries for the data referred to by the section
2694: Collective on PetscSection
2696: Input Parameters:
2697: + section - the section describing data layout
2698: - sym - the symmetry describing the affect of orientation on the access of the data
2700: Level: developer
2702: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2703: @*/
2704: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2705: {
2712: PetscObjectReference((PetscObject)sym);
2713: PetscSectionSymDestroy(&(section->sym));
2714: section->sym = sym;
2715: return(0);
2716: }
2718: /*@
2719: PetscSectionGetSym - Get the symmetries for the data referred to by the section
2721: Not collective
2723: Input Parameters:
2724: . section - the section describing data layout
2726: Output Parameters:
2727: . sym - the symmetry describing the affect of orientation on the access of the data
2729: Level: developer
2731: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2732: @*/
2733: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2734: {
2737: *sym = section->sym;
2738: return(0);
2739: }
2741: /*@
2742: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
2744: Collective on PetscSection
2746: Input Parameters:
2747: + section - the section describing data layout
2748: . field - the field number
2749: - sym - the symmetry describing the affect of orientation on the access of the data
2751: Level: developer
2753: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2754: @*/
2755: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2756: {
2761: if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
2762: PetscSectionSetSym(section->field[field],sym);
2763: return(0);
2764: }
2766: /*@
2767: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
2769: Collective on PetscSection
2771: Input Parameters:
2772: + section - the section describing data layout
2773: - field - the field number
2775: Output Parameters:
2776: . sym - the symmetry describing the affect of orientation on the access of the data
2778: Level: developer
2780: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2781: @*/
2782: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2783: {
2786: if (field < 0 || field >= section->numFields) SETERRQ2(PetscObjectComm((PetscObject)section),PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number %D (not in [0,%D)", field, section->numFields);
2787: *sym = section->field[field]->sym;
2788: return(0);
2789: }
2791: /*@C
2792: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
2794: Not collective
2796: Input Parameters:
2797: + section - the section
2798: . numPoints - the number of points
2799: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2800: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2801: context, see DMPlexGetConeOrientation()).
2803: Output Parameter:
2804: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2805: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2806: identity).
2808: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2809: .vb
2810: const PetscInt **perms;
2811: const PetscScalar **rots;
2812: PetscInt lOffset;
2814: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2815: for (i = 0, lOffset = 0; i < numPoints; i++) {
2816: PetscInt point = points[2*i], dof, sOffset;
2817: const PetscInt *perm = perms ? perms[i] : NULL;
2818: const PetscScalar *rot = rots ? rots[i] : NULL;
2820: PetscSectionGetDof(section,point,&dof);
2821: PetscSectionGetOffset(section,point,&sOffset);
2823: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
2824: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
2825: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
2826: lOffset += dof;
2827: }
2828: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2829: .ve
2831: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2832: .vb
2833: const PetscInt **perms;
2834: const PetscScalar **rots;
2835: PetscInt lOffset;
2837: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2838: for (i = 0, lOffset = 0; i < numPoints; i++) {
2839: PetscInt point = points[2*i], dof, sOffset;
2840: const PetscInt *perm = perms ? perms[i] : NULL;
2841: const PetscScalar *rot = rots ? rots[i] : NULL;
2843: PetscSectionGetDof(section,point,&dof);
2844: PetscSectionGetOffset(section,point,&sOff);
2846: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2847: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
2848: offset += dof;
2849: }
2850: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2851: .ve
2853: Level: developer
2855: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2856: @*/
2857: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2858: {
2859: PetscSectionSym sym;
2860: PetscErrorCode ierr;
2865: if (perms) *perms = NULL;
2866: if (rots) *rots = NULL;
2867: sym = section->sym;
2868: if (sym && (perms || rots)) {
2869: SymWorkLink link;
2871: if (sym->workin) {
2872: link = sym->workin;
2873: sym->workin = sym->workin->next;
2874: } else {
2875: PetscNewLog(sym,&link);
2876: }
2877: if (numPoints > link->numPoints) {
2878: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2879: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
2880: link->numPoints = numPoints;
2881: }
2882: link->next = sym->workout;
2883: sym->workout = link;
2884: PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2885: PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2886: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2887: if (perms) *perms = link->perms;
2888: if (rots) *rots = link->rots;
2889: }
2890: return(0);
2891: }
2893: /*@C
2894: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
2896: Not collective
2898: Input Parameters:
2899: + section - the section
2900: . numPoints - the number of points
2901: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2902: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2903: context, see DMPlexGetConeOrientation()).
2905: Output Parameter:
2906: + perms - The permutations for the given orientations: set to NULL at conclusion
2907: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
2909: Level: developer
2911: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2912: @*/
2913: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2914: {
2915: PetscSectionSym sym;
2919: sym = section->sym;
2920: if (sym && (perms || rots)) {
2921: SymWorkLink *p,link;
2923: for (p=&sym->workout; (link=*p); p=&link->next) {
2924: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
2925: *p = link->next;
2926: link->next = sym->workin;
2927: sym->workin = link;
2928: if (perms) *perms = NULL;
2929: if (rots) *rots = NULL;
2930: return(0);
2931: }
2932: }
2933: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
2934: }
2935: return(0);
2936: }
2938: /*@C
2939: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
2941: Not collective
2943: Input Parameters:
2944: + section - the section
2945: . field - the field of the section
2946: . numPoints - the number of points
2947: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2948: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2949: context, see DMPlexGetConeOrientation()).
2951: Output Parameter:
2952: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2953: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2954: identity).
2956: Level: developer
2958: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
2959: @*/
2960: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2961: {
2966: if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
2967: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
2968: return(0);
2969: }
2971: /*@C
2972: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
2974: Not collective
2976: Input Parameters:
2977: + section - the section
2978: . field - the field number
2979: . numPoints - the number of points
2980: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2981: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2982: context, see DMPlexGetConeOrientation()).
2984: Output Parameter:
2985: + perms - The permutations for the given orientations: set to NULL at conclusion
2986: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
2988: Level: developer
2990: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2991: @*/
2992: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2993: {
2998: if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
2999: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3000: return(0);
3001: }