Actual source code: vsectionis.c
petsc-3.10.5 2019-03-28
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)->useFieldOff = PETSC_FALSE;
60: (*s)->clObj = NULL;
61: (*s)->clSection = NULL;
62: (*s)->clPoints = NULL;
63: (*s)->clSize = 0;
64: (*s)->clPerm = NULL;
65: (*s)->clInvPerm = NULL;
66: return(0);
67: }
69: /*@
70: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
72: Collective on MPI_Comm
74: Input Parameter:
75: . section - the PetscSection
77: Output Parameter:
78: . newSection - the copy
80: Level: developer
82: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
83: @*/
84: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
85: {
86: PetscSectionSym sym;
87: IS perm;
88: PetscInt numFields, f, pStart, pEnd, p;
89: PetscErrorCode ierr;
94: PetscSectionGetNumFields(section, &numFields);
95: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
96: for (f = 0; f < numFields; ++f) {
97: const char *name = NULL;
98: PetscInt numComp = 0;
100: PetscSectionGetFieldName(section, f, &name);
101: PetscSectionSetFieldName(newSection, f, name);
102: PetscSectionGetFieldComponents(section, f, &numComp);
103: PetscSectionSetFieldComponents(newSection, f, numComp);
104: PetscSectionGetFieldSym(section, f, &sym);
105: PetscSectionSetFieldSym(newSection, f, sym);
106: }
107: PetscSectionGetChart(section, &pStart, &pEnd);
108: PetscSectionSetChart(newSection, pStart, pEnd);
109: PetscSectionGetPermutation(section, &perm);
110: PetscSectionSetPermutation(newSection, perm);
111: PetscSectionGetSym(section, &sym);
112: PetscSectionSetSym(newSection, sym);
113: for (p = pStart; p < pEnd; ++p) {
114: PetscInt dof, cdof, fcdof = 0;
116: PetscSectionGetDof(section, p, &dof);
117: PetscSectionSetDof(newSection, p, dof);
118: PetscSectionGetConstraintDof(section, p, &cdof);
119: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
120: for (f = 0; f < numFields; ++f) {
121: PetscSectionGetFieldDof(section, p, f, &dof);
122: PetscSectionSetFieldDof(newSection, p, f, dof);
123: if (cdof) {
124: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
125: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
126: }
127: }
128: }
129: PetscSectionSetUp(newSection);
130: for (p = pStart; p < pEnd; ++p) {
131: PetscInt off, cdof, fcdof = 0;
132: const PetscInt *cInd;
134: /* Must set offsets in case they do not agree with the prefix sums */
135: PetscSectionGetOffset(section, p, &off);
136: PetscSectionSetOffset(newSection, p, off);
137: PetscSectionGetConstraintDof(section, p, &cdof);
138: if (cdof) {
139: PetscSectionGetConstraintIndices(section, p, &cInd);
140: PetscSectionSetConstraintIndices(newSection, p, cInd);
141: for (f = 0; f < numFields; ++f) {
142: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
143: if (fcdof) {
144: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
145: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
146: }
147: }
148: }
149: }
150: return(0);
151: }
153: /*@
154: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
156: Collective on MPI_Comm
158: Input Parameter:
159: . section - the PetscSection
161: Output Parameter:
162: . newSection - the copy
164: Level: developer
166: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
167: @*/
168: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
169: {
175: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
176: PetscSectionCopy(section, *newSection);
177: return(0);
178: }
180: /*@
181: PetscSectionCompare - Compares two sections
183: Collective on PetscSection
185: Input Parameterss:
186: + s1 - the first PetscSection
187: - s2 - the second PetscSection
189: Output Parameter:
190: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
192: Level: developer
194: Notes:
195: Field names are disregarded.
197: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
198: @*/
199: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
200: {
201: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
202: const PetscInt *idx1, *idx2;
203: IS perm1, perm2;
204: PetscBool flg;
205: PetscMPIInt mflg;
212: flg = PETSC_FALSE;
214: MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
215: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
216: *congruent = PETSC_FALSE;
217: return(0);
218: }
220: PetscSectionGetChart(s1, &pStart, &pEnd);
221: PetscSectionGetChart(s2, &n1, &n2);
222: if (pStart != n1 || pEnd != n2) goto not_congruent;
224: PetscSectionGetPermutation(s1, &perm1);
225: PetscSectionGetPermutation(s2, &perm2);
226: if (perm1 && perm2) {
227: ISEqual(perm1, perm2, congruent);
228: if (!(*congruent)) goto not_congruent;
229: } else if (perm1 != perm2) goto not_congruent;
231: for (p = pStart; p < pEnd; ++p) {
232: PetscSectionGetOffset(s1, p, &n1);
233: PetscSectionGetOffset(s2, p, &n2);
234: if (n1 != n2) goto not_congruent;
236: PetscSectionGetDof(s1, p, &n1);
237: PetscSectionGetDof(s2, p, &n2);
238: if (n1 != n2) goto not_congruent;
240: PetscSectionGetConstraintDof(s1, p, &ncdof);
241: PetscSectionGetConstraintDof(s2, p, &n2);
242: if (ncdof != n2) goto not_congruent;
244: PetscSectionGetConstraintIndices(s1, p, &idx1);
245: PetscSectionGetConstraintIndices(s2, p, &idx2);
246: PetscMemcmp(idx1, idx2, ncdof*sizeof(PetscInt), congruent);
247: if (!(*congruent)) goto not_congruent;
248: }
250: PetscSectionGetNumFields(s1, &nfields);
251: PetscSectionGetNumFields(s2, &n2);
252: if (nfields != n2) goto not_congruent;
254: for (f = 0; f < nfields; ++f) {
255: PetscSectionGetFieldComponents(s1, f, &n1);
256: PetscSectionGetFieldComponents(s2, f, &n2);
257: if (n1 != n2) goto not_congruent;
259: for (p = pStart; p < pEnd; ++p) {
260: PetscSectionGetFieldOffset(s1, p, f, &n1);
261: PetscSectionGetFieldOffset(s2, p, f, &n2);
262: if (n1 != n2) goto not_congruent;
264: PetscSectionGetFieldDof(s1, p, f, &n1);
265: PetscSectionGetFieldDof(s2, p, f, &n2);
266: if (n1 != n2) goto not_congruent;
268: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
269: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
270: if (nfcdof != n2) goto not_congruent;
272: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
273: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
274: PetscMemcmp(idx1, idx2, nfcdof*sizeof(PetscInt), congruent);
275: if (!(*congruent)) goto not_congruent;
276: }
277: }
279: flg = PETSC_TRUE;
280: not_congruent:
281: MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
282: return(0);
283: }
285: /*@
286: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
288: Not collective
290: Input Parameter:
291: . s - the PetscSection
293: Output Parameter:
294: . numFields - the number of fields defined, or 0 if none were defined
296: Level: intermediate
298: .seealso: PetscSectionSetNumFields()
299: @*/
300: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
301: {
305: *numFields = s->numFields;
306: return(0);
307: }
309: /*@
310: PetscSectionSetNumFields - Sets the number of fields.
312: Not collective
314: Input Parameters:
315: + s - the PetscSection
316: - numFields - the number of fields
318: Level: intermediate
320: .seealso: PetscSectionGetNumFields()
321: @*/
322: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
323: {
324: PetscInt f;
329: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %d must be positive", numFields);
330: PetscSectionReset(s);
332: s->numFields = numFields;
333: PetscMalloc1(s->numFields, &s->numFieldComponents);
334: PetscMalloc1(s->numFields, &s->fieldNames);
335: PetscMalloc1(s->numFields, &s->field);
336: for (f = 0; f < s->numFields; ++f) {
337: char name[64];
339: s->numFieldComponents[f] = 1;
341: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
342: PetscSNPrintf(name, 64, "Field_%D", f);
343: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
344: }
345: return(0);
346: }
348: /*@C
349: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
351: Not Collective
353: Input Parameters:
354: + s - the PetscSection
355: - field - the field number
357: Output Parameter:
358: . fieldName - the field name
360: Level: developer
362: .seealso: PetscSectionSetFieldName()
363: @*/
364: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
365: {
369: 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);
370: *fieldName = s->fieldNames[field];
371: return(0);
372: }
374: /*@C
375: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
377: Not Collective
379: Input Parameters:
380: + s - the PetscSection
381: . field - the field number
382: - fieldName - the field name
384: Level: developer
386: .seealso: PetscSectionGetFieldName()
387: @*/
388: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
389: {
395: 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);
396: PetscFree(s->fieldNames[field]);
397: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
398: return(0);
399: }
401: /*@
402: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
404: Not collective
406: Input Parameters:
407: + s - the PetscSection
408: - field - the field number
410: Output Parameter:
411: . numComp - the number of field components
413: Level: intermediate
415: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
416: @*/
417: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
418: {
422: 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);
423: *numComp = s->numFieldComponents[field];
424: return(0);
425: }
427: /*@
428: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
430: Not collective
432: Input Parameters:
433: + s - the PetscSection
434: . field - the field number
435: - numComp - the number of field components
437: Level: intermediate
439: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
440: @*/
441: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
442: {
445: 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);
446: s->numFieldComponents[field] = numComp;
447: return(0);
448: }
450: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
451: {
455: if (!s->bc) {
456: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
457: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
458: }
459: return(0);
460: }
462: /*@
463: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
465: Not collective
467: Input Parameter:
468: . s - the PetscSection
470: Output Parameters:
471: + pStart - the first point
472: - pEnd - one past the last point
474: Level: intermediate
476: .seealso: PetscSectionSetChart(), PetscSectionCreate()
477: @*/
478: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
479: {
482: if (pStart) *pStart = s->pStart;
483: if (pEnd) *pEnd = s->pEnd;
484: return(0);
485: }
487: /*@
488: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
490: Not collective
492: Input Parameters:
493: + s - the PetscSection
494: . pStart - the first point
495: - pEnd - one past the last point
497: Level: intermediate
499: .seealso: PetscSectionGetChart(), PetscSectionCreate()
500: @*/
501: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
502: {
503: PetscInt f;
508: /* Cannot Reset() because it destroys field information */
509: s->setup = PETSC_FALSE;
510: PetscSectionDestroy(&s->bc);
511: PetscFree(s->bcIndices);
512: PetscFree2(s->atlasDof, s->atlasOff);
514: s->pStart = pStart;
515: s->pEnd = pEnd;
516: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
517: PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
518: for (f = 0; f < s->numFields; ++f) {
519: PetscSectionSetChart(s->field[f], pStart, pEnd);
520: }
521: return(0);
522: }
524: /*@
525: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
527: Not collective
529: Input Parameter:
530: . s - the PetscSection
532: Output Parameters:
533: . perm - The permutation as an IS
535: Level: intermediate
537: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
538: @*/
539: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
540: {
544: return(0);
545: }
547: /*@
548: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
550: Not collective
552: Input Parameters:
553: + s - the PetscSection
554: - perm - the permutation of points
556: Level: intermediate
558: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
559: @*/
560: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
561: {
567: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
568: if (s->perm != perm) {
569: ISDestroy(&s->perm);
570: if (perm) {
571: s->perm = perm;
572: PetscObjectReference((PetscObject) s->perm);
573: }
574: }
575: return(0);
576: }
578: /*@
579: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
581: Not collective
583: Input Parameters:
584: + s - the PetscSection
585: - point - the point
587: Output Parameter:
588: . numDof - the number of dof
590: Level: intermediate
592: .seealso: PetscSectionSetDof(), PetscSectionCreate()
593: @*/
594: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
595: {
599: #if defined(PETSC_USE_DEBUG)
600: 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);
601: #endif
602: *numDof = s->atlasDof[point - s->pStart];
603: return(0);
604: }
606: /*@
607: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
609: Not collective
611: Input Parameters:
612: + s - the PetscSection
613: . point - the point
614: - numDof - the number of dof
616: Level: intermediate
618: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
619: @*/
620: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
621: {
624: 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);
625: s->atlasDof[point - s->pStart] = numDof;
626: return(0);
627: }
629: /*@
630: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
632: Not collective
634: Input Parameters:
635: + s - the PetscSection
636: . point - the point
637: - numDof - the number of additional dof
639: Level: intermediate
641: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
642: @*/
643: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
644: {
647: 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);
648: s->atlasDof[point - s->pStart] += numDof;
649: return(0);
650: }
652: /*@
653: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
655: Not collective
657: Input Parameters:
658: + s - the PetscSection
659: . point - the point
660: - field - the field
662: Output Parameter:
663: . numDof - the number of dof
665: Level: intermediate
667: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
668: @*/
669: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
670: {
675: 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);
676: PetscSectionGetDof(s->field[field], point, numDof);
677: return(0);
678: }
680: /*@
681: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
683: Not collective
685: Input Parameters:
686: + s - the PetscSection
687: . point - the point
688: . field - the field
689: - numDof - the number of dof
691: Level: intermediate
693: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
694: @*/
695: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
696: {
701: 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);
702: PetscSectionSetDof(s->field[field], point, numDof);
703: return(0);
704: }
706: /*@
707: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
709: Not collective
711: Input Parameters:
712: + s - the PetscSection
713: . point - the point
714: . field - the field
715: - numDof - the number of dof
717: Level: intermediate
719: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
720: @*/
721: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
722: {
727: 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);
728: PetscSectionAddDof(s->field[field], point, numDof);
729: return(0);
730: }
732: /*@
733: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
735: Not collective
737: Input Parameters:
738: + s - the PetscSection
739: - point - the point
741: Output Parameter:
742: . numDof - the number of dof which are fixed by constraints
744: Level: intermediate
746: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
747: @*/
748: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
749: {
755: if (s->bc) {
756: PetscSectionGetDof(s->bc, point, numDof);
757: } else *numDof = 0;
758: return(0);
759: }
761: /*@
762: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
764: Not collective
766: Input Parameters:
767: + s - the PetscSection
768: . point - the point
769: - numDof - the number of dof which are fixed by constraints
771: Level: intermediate
773: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
774: @*/
775: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
776: {
781: if (numDof) {
782: PetscSectionCheckConstraints_Static(s);
783: PetscSectionSetDof(s->bc, point, numDof);
784: }
785: return(0);
786: }
788: /*@
789: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
791: Not collective
793: Input Parameters:
794: + s - the PetscSection
795: . point - the point
796: - numDof - the number of additional dof which are fixed by constraints
798: Level: intermediate
800: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
801: @*/
802: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
803: {
808: if (numDof) {
809: PetscSectionCheckConstraints_Static(s);
810: PetscSectionAddDof(s->bc, point, numDof);
811: }
812: return(0);
813: }
815: /*@
816: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
818: Not collective
820: Input Parameters:
821: + s - the PetscSection
822: . point - the point
823: - field - the field
825: Output Parameter:
826: . numDof - the number of dof which are fixed by constraints
828: Level: intermediate
830: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
831: @*/
832: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
833: {
839: 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);
840: PetscSectionGetConstraintDof(s->field[field], point, numDof);
841: return(0);
842: }
844: /*@
845: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
847: Not collective
849: Input Parameters:
850: + s - the PetscSection
851: . point - the point
852: . field - the field
853: - numDof - the number of dof which are fixed by constraints
855: Level: intermediate
857: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
858: @*/
859: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
860: {
865: 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);
866: PetscSectionSetConstraintDof(s->field[field], point, numDof);
867: return(0);
868: }
870: /*@
871: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
873: Not collective
875: Input Parameters:
876: + s - the PetscSection
877: . point - the point
878: . field - the field
879: - numDof - the number of additional dof which are fixed by constraints
881: Level: intermediate
883: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
884: @*/
885: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
886: {
891: 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);
892: PetscSectionAddConstraintDof(s->field[field], point, numDof);
893: return(0);
894: }
896: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
897: {
902: if (s->bc) {
903: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
905: PetscSectionSetUp(s->bc);
906: PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
907: }
908: return(0);
909: }
911: /*@
912: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
914: Not collective
916: Input Parameter:
917: . s - the PetscSection
919: Level: intermediate
921: .seealso: PetscSectionCreate()
922: @*/
923: PetscErrorCode PetscSectionSetUp(PetscSection s)
924: {
925: const PetscInt *pind = NULL;
926: PetscInt offset = 0, p, f;
927: PetscErrorCode ierr;
931: if (s->setup) return(0);
932: s->setup = PETSC_TRUE;
933: if (s->perm) {ISGetIndices(s->perm, &pind);}
934: for (p = 0; p < s->pEnd - s->pStart; ++p) {
935: const PetscInt q = pind ? pind[p] : p;
937: s->atlasOff[q] = offset;
938: offset += s->atlasDof[q];
939: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
940: }
941: PetscSectionSetUpBC(s);
942: /* Assume that all fields have the same chart */
943: for (p = 0; p < s->pEnd - s->pStart; ++p) {
944: const PetscInt q = pind ? pind[p] : p;
945: PetscInt off = s->atlasOff[q];
947: for (f = 0; f < s->numFields; ++f) {
948: PetscSection sf = s->field[f];
950: sf->atlasOff[q] = off;
951: off += sf->atlasDof[q];
952: }
953: }
954: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
955: for (f = 0; f < s->numFields; ++f) {
956: PetscSectionSetUpBC(s->field[f]);
957: }
958: return(0);
959: }
961: /*@
962: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
964: Not collective
966: Input Parameters:
967: . s - the PetscSection
969: Output Parameter:
970: . maxDof - the maximum dof
972: Level: intermediate
974: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
975: @*/
976: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
977: {
981: *maxDof = s->maxDof;
982: return(0);
983: }
985: /*@
986: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
988: Not collective
990: Input Parameters:
991: + s - the PetscSection
992: - size - the allocated size
994: Output Parameter:
995: . size - the size of an array which can hold all the dofs
997: Level: intermediate
999: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1000: @*/
1001: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1002: {
1003: PetscInt p, n = 0;
1008: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1009: *size = n;
1010: return(0);
1011: }
1013: /*@
1014: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1016: Not collective
1018: Input Parameters:
1019: + s - the PetscSection
1020: - point - the point
1022: Output Parameter:
1023: . size - the size of an array which can hold all unconstrained dofs
1025: Level: intermediate
1027: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1028: @*/
1029: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1030: {
1031: PetscInt p, n = 0;
1036: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1037: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1038: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1039: }
1040: *size = n;
1041: return(0);
1042: }
1044: /*@
1045: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1046: the local section and an SF describing the section point overlap.
1048: Input Parameters:
1049: + s - The PetscSection for the local field layout
1050: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1051: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1052: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1054: Output Parameter:
1055: . gsection - The PetscSection for the global field layout
1057: Note: This gives negative sizes and offsets to points not owned by this process
1059: Level: developer
1061: .seealso: PetscSectionCreate()
1062: @*/
1063: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1064: {
1065: PetscSection gs;
1066: const PetscInt *pind = NULL;
1067: PetscInt *recv = NULL, *neg = NULL;
1068: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
1069: PetscErrorCode ierr;
1077: PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1078: PetscSectionGetChart(s, &pStart, &pEnd);
1079: PetscSectionSetChart(gs, pStart, pEnd);
1080: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1081: nlocal = nroots; /* The local/leaf space matches global/root space */
1082: /* Must allocate for all points visible to SF, which may be more than this section */
1083: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1084: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
1085: PetscMalloc2(nroots,&neg,nlocal,&recv);
1086: PetscMemzero(neg,nroots*sizeof(PetscInt));
1087: }
1088: /* Mark all local points with negative dof */
1089: for (p = pStart; p < pEnd; ++p) {
1090: PetscSectionGetDof(s, p, &dof);
1091: PetscSectionSetDof(gs, p, dof);
1092: PetscSectionGetConstraintDof(s, p, &cdof);
1093: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1094: if (neg) neg[p] = -(dof+1);
1095: }
1096: PetscSectionSetUpBC(gs);
1097: 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]));}
1098: if (nroots >= 0) {
1099: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1100: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1101: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1102: for (p = pStart; p < pEnd; ++p) {
1103: if (recv[p] < 0) {
1104: gs->atlasDof[p-pStart] = recv[p];
1105: PetscSectionGetDof(s, p, &dof);
1106: 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);
1107: }
1108: }
1109: }
1110: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1111: if (s->perm) {ISGetIndices(s->perm, &pind);}
1112: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1113: const PetscInt q = pind ? pind[p] : p;
1115: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1116: gs->atlasOff[q] = off;
1117: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1118: }
1119: if (!localOffsets) {
1120: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1121: globalOff -= off;
1122: }
1123: for (p = pStart, off = 0; p < pEnd; ++p) {
1124: gs->atlasOff[p-pStart] += globalOff;
1125: if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1126: }
1127: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1128: /* Put in negative offsets for ghost points */
1129: if (nroots >= 0) {
1130: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1131: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1132: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1133: for (p = pStart; p < pEnd; ++p) {
1134: if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1135: }
1136: }
1137: PetscFree2(neg,recv);
1138: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1139: *gsection = gs;
1140: return(0);
1141: }
1143: /*@
1144: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1145: the local section and an SF describing the section point overlap.
1147: Input Parameters:
1148: + s - The PetscSection for the local field layout
1149: . sf - The SF describing parallel layout of the section points
1150: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1151: . numExcludes - The number of exclusion ranges
1152: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1154: Output Parameter:
1155: . gsection - The PetscSection for the global field layout
1157: Note: This gives negative sizes and offsets to points not owned by this process
1159: Level: developer
1161: .seealso: PetscSectionCreate()
1162: @*/
1163: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1164: {
1165: const PetscInt *pind = NULL;
1166: PetscInt *neg = NULL, *tmpOff = NULL;
1167: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1168: PetscErrorCode ierr;
1174: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1175: PetscSectionGetChart(s, &pStart, &pEnd);
1176: PetscSectionSetChart(*gsection, pStart, pEnd);
1177: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1178: if (nroots >= 0) {
1179: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1180: PetscCalloc1(nroots, &neg);
1181: if (nroots > pEnd-pStart) {
1182: PetscCalloc1(nroots, &tmpOff);
1183: } else {
1184: tmpOff = &(*gsection)->atlasDof[-pStart];
1185: }
1186: }
1187: /* Mark ghost points with negative dof */
1188: for (p = pStart; p < pEnd; ++p) {
1189: for (e = 0; e < numExcludes; ++e) {
1190: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1191: PetscSectionSetDof(*gsection, p, 0);
1192: break;
1193: }
1194: }
1195: if (e < numExcludes) continue;
1196: PetscSectionGetDof(s, p, &dof);
1197: PetscSectionSetDof(*gsection, p, dof);
1198: PetscSectionGetConstraintDof(s, p, &cdof);
1199: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1200: if (neg) neg[p] = -(dof+1);
1201: }
1202: PetscSectionSetUpBC(*gsection);
1203: if (nroots >= 0) {
1204: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1205: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1206: if (nroots > pEnd - pStart) {
1207: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1208: }
1209: }
1210: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1211: if (s->perm) {ISGetIndices(s->perm, &pind);}
1212: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1213: const PetscInt q = pind ? pind[p] : p;
1215: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1216: (*gsection)->atlasOff[q] = off;
1217: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1218: }
1219: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1220: globalOff -= off;
1221: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1222: (*gsection)->atlasOff[p] += globalOff;
1223: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1224: }
1225: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1226: /* Put in negative offsets for ghost points */
1227: if (nroots >= 0) {
1228: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1229: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1230: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1231: if (nroots > pEnd - pStart) {
1232: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1233: }
1234: }
1235: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1236: PetscFree(neg);
1237: return(0);
1238: }
1240: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1241: {
1242: PetscInt pStart, pEnd, p, localSize = 0;
1246: PetscSectionGetChart(s, &pStart, &pEnd);
1247: for (p = pStart; p < pEnd; ++p) {
1248: PetscInt dof;
1250: PetscSectionGetDof(s, p, &dof);
1251: if (dof > 0) ++localSize;
1252: }
1253: PetscLayoutCreate(comm, layout);
1254: PetscLayoutSetLocalSize(*layout, localSize);
1255: PetscLayoutSetBlockSize(*layout, 1);
1256: PetscLayoutSetUp(*layout);
1257: return(0);
1258: }
1260: /*@
1261: PetscSectionGetValueLayout - Get the PetscLayout associated with a section, usually the default global section.
1263: Input Parameters:
1264: + comm - The MPI_Comm
1265: - s - The PetscSection
1267: Output Parameter:
1268: . layout - The layout for the section
1270: Level: developer
1272: .seealso: PetscSectionCreate()
1273: @*/
1274: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1275: {
1276: PetscInt pStart, pEnd, p, localSize = 0;
1282: PetscSectionGetChart(s, &pStart, &pEnd);
1283: for (p = pStart; p < pEnd; ++p) {
1284: PetscInt dof,cdof;
1286: PetscSectionGetDof(s, p, &dof);
1287: PetscSectionGetConstraintDof(s, p, &cdof);
1288: if (dof-cdof > 0) localSize += dof-cdof;
1289: }
1290: PetscLayoutCreate(comm, layout);
1291: PetscLayoutSetLocalSize(*layout, localSize);
1292: PetscLayoutSetBlockSize(*layout, 1);
1293: PetscLayoutSetUp(*layout);
1294: return(0);
1295: }
1297: /*@
1298: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1300: Not collective
1302: Input Parameters:
1303: + s - the PetscSection
1304: - point - the point
1306: Output Parameter:
1307: . offset - the offset
1309: Level: intermediate
1311: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1312: @*/
1313: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1314: {
1318: #if defined(PETSC_USE_DEBUG)
1319: 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);
1320: #endif
1321: *offset = s->atlasOff[point - s->pStart];
1322: return(0);
1323: }
1325: /*@
1326: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1328: Not collective
1330: Input Parameters:
1331: + s - the PetscSection
1332: . point - the point
1333: - offset - the offset
1335: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1337: Level: intermediate
1339: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1340: @*/
1341: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1342: {
1345: 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);
1346: s->atlasOff[point - s->pStart] = offset;
1347: return(0);
1348: }
1350: /*@
1351: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1353: Not collective
1355: Input Parameters:
1356: + s - the PetscSection
1357: . point - the point
1358: - field - the field
1360: Output Parameter:
1361: . offset - the offset
1363: Level: intermediate
1365: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1366: @*/
1367: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1368: {
1374: 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);
1375: PetscSectionGetOffset(s->field[field], point, offset);
1376: return(0);
1377: }
1379: /*@
1380: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1382: Not collective
1384: Input Parameters:
1385: + s - the PetscSection
1386: . point - the point
1387: . field - the field
1388: - offset - the offset
1390: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1392: Level: intermediate
1394: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1395: @*/
1396: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1397: {
1402: 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);
1403: PetscSectionSetOffset(s->field[field], point, offset);
1404: return(0);
1405: }
1407: /* This gives the offset on a point of the field, ignoring constraints */
1408: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1409: {
1410: PetscInt off, foff;
1416: 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);
1417: PetscSectionGetOffset(s, point, &off);
1418: PetscSectionGetOffset(s->field[field], point, &foff);
1419: *offset = foff - off;
1420: return(0);
1421: }
1423: /*@
1424: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1426: Not collective
1428: Input Parameter:
1429: . s - the PetscSection
1431: Output Parameters:
1432: + start - the minimum offset
1433: - end - one more than the maximum offset
1435: Level: intermediate
1437: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1438: @*/
1439: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1440: {
1441: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1446: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1447: PetscSectionGetChart(s, &pStart, &pEnd);
1448: for (p = 0; p < pEnd-pStart; ++p) {
1449: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1451: if (off >= 0) {
1452: os = PetscMin(os, off);
1453: oe = PetscMax(oe, off+dof);
1454: }
1455: }
1456: if (start) *start = os;
1457: if (end) *end = oe;
1458: return(0);
1459: }
1461: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, const PetscInt fields[], PetscSection *subs)
1462: {
1463: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1467: if (!numFields) return(0);
1471: PetscSectionGetNumFields(s, &nF);
1472: if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1473: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1474: PetscSectionSetNumFields(*subs, numFields);
1475: for (f = 0; f < numFields; ++f) {
1476: const char *name = NULL;
1477: PetscInt numComp = 0;
1479: PetscSectionGetFieldName(s, fields[f], &name);
1480: PetscSectionSetFieldName(*subs, f, name);
1481: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1482: PetscSectionSetFieldComponents(*subs, f, numComp);
1483: }
1484: PetscSectionGetChart(s, &pStart, &pEnd);
1485: PetscSectionSetChart(*subs, pStart, pEnd);
1486: for (p = pStart; p < pEnd; ++p) {
1487: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1489: for (f = 0; f < numFields; ++f) {
1490: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1491: PetscSectionSetFieldDof(*subs, p, f, fdof);
1492: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1493: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1494: dof += fdof;
1495: cdof += cfdof;
1496: }
1497: PetscSectionSetDof(*subs, p, dof);
1498: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1499: maxCdof = PetscMax(cdof, maxCdof);
1500: }
1501: PetscSectionSetUp(*subs);
1502: if (maxCdof) {
1503: PetscInt *indices;
1505: PetscMalloc1(maxCdof, &indices);
1506: for (p = pStart; p < pEnd; ++p) {
1507: PetscInt cdof;
1509: PetscSectionGetConstraintDof(*subs, p, &cdof);
1510: if (cdof) {
1511: const PetscInt *oldIndices = NULL;
1512: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1514: for (f = 0; f < numFields; ++f) {
1515: PetscInt oldFoff = 0, oldf;
1517: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1518: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1519: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1520: /* This can be sped up if we assume sorted fields */
1521: for (oldf = 0; oldf < fields[f]; ++oldf) {
1522: PetscInt oldfdof = 0;
1523: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1524: oldFoff += oldfdof;
1525: }
1526: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1527: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1528: numConst += cfdof;
1529: fOff += fdof;
1530: }
1531: PetscSectionSetConstraintIndices(*subs, p, indices);
1532: }
1533: }
1534: PetscFree(indices);
1535: }
1536: return(0);
1537: }
1539: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1540: {
1541: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1545: if (!len) return(0);
1546: for (i = 0; i < len; ++i) {
1547: PetscInt nf, pStarti, pEndi;
1549: PetscSectionGetNumFields(s[i], &nf);
1550: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1551: pStart = PetscMin(pStart, pStarti);
1552: pEnd = PetscMax(pEnd, pEndi);
1553: Nf += nf;
1554: }
1555: PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1556: PetscSectionSetNumFields(*supers, Nf);
1557: for (i = 0, f = 0; i < len; ++i) {
1558: PetscInt nf, fi;
1560: PetscSectionGetNumFields(s[i], &nf);
1561: for (fi = 0; fi < nf; ++fi, ++f) {
1562: const char *name = NULL;
1563: PetscInt numComp = 0;
1565: PetscSectionGetFieldName(s[i], fi, &name);
1566: PetscSectionSetFieldName(*supers, f, name);
1567: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1568: PetscSectionSetFieldComponents(*supers, f, numComp);
1569: }
1570: }
1571: PetscSectionSetChart(*supers, pStart, pEnd);
1572: for (p = pStart; p < pEnd; ++p) {
1573: PetscInt dof = 0, cdof = 0;
1575: for (i = 0, f = 0; i < len; ++i) {
1576: PetscInt nf, fi, pStarti, pEndi;
1577: PetscInt fdof = 0, cfdof = 0;
1579: PetscSectionGetNumFields(s[i], &nf);
1580: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1581: if ((p < pStarti) || (p >= pEndi)) continue;
1582: for (fi = 0; fi < nf; ++fi, ++f) {
1583: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1584: PetscSectionAddFieldDof(*supers, p, f, fdof);
1585: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1586: if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1587: dof += fdof;
1588: cdof += cfdof;
1589: }
1590: }
1591: PetscSectionSetDof(*supers, p, dof);
1592: if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1593: maxCdof = PetscMax(cdof, maxCdof);
1594: }
1595: PetscSectionSetUp(*supers);
1596: if (maxCdof) {
1597: PetscInt *indices;
1599: PetscMalloc1(maxCdof, &indices);
1600: for (p = pStart; p < pEnd; ++p) {
1601: PetscInt cdof;
1603: PetscSectionGetConstraintDof(*supers, p, &cdof);
1604: if (cdof) {
1605: PetscInt dof, numConst = 0, fOff = 0;
1607: for (i = 0, f = 0; i < len; ++i) {
1608: const PetscInt *oldIndices = NULL;
1609: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1611: PetscSectionGetNumFields(s[i], &nf);
1612: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1613: if ((p < pStarti) || (p >= pEndi)) continue;
1614: for (fi = 0; fi < nf; ++fi, ++f) {
1615: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1616: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1617: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1618: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1619: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1620: numConst += cfdof;
1621: }
1622: PetscSectionGetDof(s[i], p, &dof);
1623: fOff += dof;
1624: }
1625: PetscSectionSetConstraintIndices(*supers, p, indices);
1626: }
1627: }
1628: PetscFree(indices);
1629: }
1630: return(0);
1631: }
1633: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1634: {
1635: const PetscInt *points = NULL, *indices = NULL;
1636: PetscInt numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;
1643: PetscSectionGetNumFields(s, &numFields);
1644: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1645: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1646: for (f = 0; f < numFields; ++f) {
1647: const char *name = NULL;
1648: PetscInt numComp = 0;
1650: PetscSectionGetFieldName(s, f, &name);
1651: PetscSectionSetFieldName(*subs, f, name);
1652: PetscSectionGetFieldComponents(s, f, &numComp);
1653: PetscSectionSetFieldComponents(*subs, f, numComp);
1654: }
1655: /* For right now, we do not try to squeeze the subchart */
1656: if (subpointMap) {
1657: ISGetSize(subpointMap, &numSubpoints);
1658: ISGetIndices(subpointMap, &points);
1659: }
1660: PetscSectionGetChart(s, &pStart, &pEnd);
1661: PetscSectionSetChart(*subs, 0, numSubpoints);
1662: for (p = pStart; p < pEnd; ++p) {
1663: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1665: PetscFindInt(p, numSubpoints, points, &subp);
1666: if (subp < 0) continue;
1667: for (f = 0; f < numFields; ++f) {
1668: PetscSectionGetFieldDof(s, p, f, &fdof);
1669: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1670: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1671: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1672: }
1673: PetscSectionGetDof(s, p, &dof);
1674: PetscSectionSetDof(*subs, subp, dof);
1675: PetscSectionGetConstraintDof(s, p, &cdof);
1676: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1677: }
1678: PetscSectionSetUp(*subs);
1679: /* Change offsets to original offsets */
1680: for (p = pStart; p < pEnd; ++p) {
1681: PetscInt off, foff = 0;
1683: PetscFindInt(p, numSubpoints, points, &subp);
1684: if (subp < 0) continue;
1685: for (f = 0; f < numFields; ++f) {
1686: PetscSectionGetFieldOffset(s, p, f, &foff);
1687: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1688: }
1689: PetscSectionGetOffset(s, p, &off);
1690: PetscSectionSetOffset(*subs, subp, off);
1691: }
1692: /* Copy constraint indices */
1693: for (subp = 0; subp < numSubpoints; ++subp) {
1694: PetscInt cdof;
1696: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1697: if (cdof) {
1698: for (f = 0; f < numFields; ++f) {
1699: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1700: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1701: }
1702: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1703: PetscSectionSetConstraintIndices(*subs, subp, indices);
1704: }
1705: }
1706: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1707: return(0);
1708: }
1710: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1711: {
1712: PetscInt p;
1713: PetscMPIInt rank;
1717: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1718: PetscViewerASCIIPushSynchronized(viewer);
1719: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1720: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1721: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1722: PetscInt b;
1724: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1725: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1726: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1727: }
1728: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1729: } else {
1730: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1731: }
1732: }
1733: PetscViewerFlush(viewer);
1734: PetscViewerASCIIPopSynchronized(viewer);
1735: if (s->sym) {
1736: PetscViewerASCIIPushTab(viewer);
1737: PetscSectionSymView(s->sym,viewer);
1738: PetscViewerASCIIPopTab(viewer);
1739: }
1740: return(0);
1741: }
1743: /*@C
1744: PetscSectionView - Views a PetscSection
1746: Collective on PetscSection
1748: Input Parameters:
1749: + s - the PetscSection object to view
1750: - v - the viewer
1752: Level: developer
1754: .seealso PetscSectionCreate(), PetscSectionDestroy()
1755: @*/
1756: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1757: {
1758: PetscBool isascii;
1759: PetscInt f;
1764: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1766: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1767: if (isascii) {
1768: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1769: if (s->numFields) {
1770: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1771: for (f = 0; f < s->numFields; ++f) {
1772: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
1773: PetscSectionView_ASCII(s->field[f], viewer);
1774: }
1775: } else {
1776: PetscSectionView_ASCII(s, viewer);
1777: }
1778: }
1779: return(0);
1780: }
1782: /*@
1783: PetscSectionReset - Frees all section data.
1785: Not collective
1787: Input Parameters:
1788: . s - the PetscSection
1790: Level: developer
1792: .seealso: PetscSection, PetscSectionCreate()
1793: @*/
1794: PetscErrorCode PetscSectionReset(PetscSection s)
1795: {
1796: PetscInt f;
1801: PetscFree(s->numFieldComponents);
1802: for (f = 0; f < s->numFields; ++f) {
1803: PetscSectionDestroy(&s->field[f]);
1804: PetscFree(s->fieldNames[f]);
1805: }
1806: PetscFree(s->fieldNames);
1807: PetscFree(s->field);
1808: PetscSectionDestroy(&s->bc);
1809: PetscFree(s->bcIndices);
1810: PetscFree2(s->atlasDof, s->atlasOff);
1811: PetscSectionDestroy(&s->clSection);
1812: ISDestroy(&s->clPoints);
1813: ISDestroy(&s->perm);
1814: PetscFree(s->clPerm);
1815: PetscFree(s->clInvPerm);
1816: PetscSectionSymDestroy(&s->sym);
1818: s->pStart = -1;
1819: s->pEnd = -1;
1820: s->maxDof = 0;
1821: s->setup = PETSC_FALSE;
1822: s->numFields = 0;
1823: s->clObj = NULL;
1824: return(0);
1825: }
1827: /*@
1828: PetscSectionDestroy - Frees a section object and frees its range if that exists.
1830: Not collective
1832: Input Parameters:
1833: . s - the PetscSection
1835: Level: developer
1837: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1838: recommended they not be used in user codes unless you really gain something in their use.
1840: .seealso: PetscSection, PetscSectionCreate()
1841: @*/
1842: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1843: {
1847: if (!*s) return(0);
1849: if (--((PetscObject)(*s))->refct > 0) {
1850: *s = NULL;
1851: return(0);
1852: }
1853: PetscSectionReset(*s);
1854: PetscHeaderDestroy(s);
1855: return(0);
1856: }
1858: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1859: {
1860: const PetscInt p = point - s->pStart;
1864: *values = &baseArray[s->atlasOff[p]];
1865: return(0);
1866: }
1868: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1869: {
1870: PetscInt *array;
1871: const PetscInt p = point - s->pStart;
1872: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1873: PetscInt cDim = 0;
1878: PetscSectionGetConstraintDof(s, p, &cDim);
1879: array = &baseArray[s->atlasOff[p]];
1880: if (!cDim) {
1881: if (orientation >= 0) {
1882: const PetscInt dim = s->atlasDof[p];
1883: PetscInt i;
1885: if (mode == INSERT_VALUES) {
1886: for (i = 0; i < dim; ++i) array[i] = values[i];
1887: } else {
1888: for (i = 0; i < dim; ++i) array[i] += values[i];
1889: }
1890: } else {
1891: PetscInt offset = 0;
1892: PetscInt j = -1, field, i;
1894: for (field = 0; field < s->numFields; ++field) {
1895: const PetscInt dim = s->field[field]->atlasDof[p];
1897: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1898: offset += dim;
1899: }
1900: }
1901: } else {
1902: if (orientation >= 0) {
1903: const PetscInt dim = s->atlasDof[p];
1904: PetscInt cInd = 0, i;
1905: const PetscInt *cDof;
1907: PetscSectionGetConstraintIndices(s, point, &cDof);
1908: if (mode == INSERT_VALUES) {
1909: for (i = 0; i < dim; ++i) {
1910: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1911: array[i] = values[i];
1912: }
1913: } else {
1914: for (i = 0; i < dim; ++i) {
1915: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1916: array[i] += values[i];
1917: }
1918: }
1919: } else {
1920: const PetscInt *cDof;
1921: PetscInt offset = 0;
1922: PetscInt cOffset = 0;
1923: PetscInt j = 0, field;
1925: PetscSectionGetConstraintIndices(s, point, &cDof);
1926: for (field = 0; field < s->numFields; ++field) {
1927: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
1928: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1929: const PetscInt sDim = dim - tDim;
1930: PetscInt cInd = 0, i,k;
1932: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1933: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1934: array[j] = values[k];
1935: }
1936: offset += dim;
1937: cOffset += dim - tDim;
1938: }
1939: }
1940: }
1941: return(0);
1942: }
1944: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1945: {
1949: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1950: return(0);
1951: }
1953: /*@C
1954: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
1956: Input Parameters:
1957: + s - The PetscSection
1958: - point - The point
1960: Output Parameter:
1961: . indices - The constrained dofs
1963: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
1965: Level: advanced
1967: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1968: @*/
1969: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1970: {
1975: if (s->bc) {
1976: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1977: } else *indices = NULL;
1978: return(0);
1979: }
1981: /*@C
1982: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
1984: Input Parameters:
1985: + s - The PetscSection
1986: . point - The point
1987: - indices - The constrained dofs
1989: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
1991: Level: advanced
1993: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1994: @*/
1995: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1996: {
2001: if (s->bc) {
2002: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2003: }
2004: return(0);
2005: }
2007: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2008: {
2013: 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);
2014: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2015: return(0);
2016: }
2018: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2019: {
2024: 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);
2025: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2026: return(0);
2027: }
2029: /*@
2030: PetscSectionPermute - Reorder the section according to the input point permutation
2032: Collective on PetscSection
2034: Input Parameter:
2035: + section - The PetscSection object
2036: - perm - The point permutation, old point p becomes new point perm[p]
2038: Output Parameter:
2039: . sectionNew - The permuted PetscSection
2041: Level: intermediate
2043: .keywords: mesh
2044: .seealso: MatPermute()
2045: @*/
2046: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2047: {
2048: PetscSection s = section, sNew;
2049: const PetscInt *perm;
2050: PetscInt numFields, f, numPoints, pStart, pEnd, p;
2051: PetscErrorCode ierr;
2057: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2058: PetscSectionGetNumFields(s, &numFields);
2059: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2060: for (f = 0; f < numFields; ++f) {
2061: const char *name;
2062: PetscInt numComp;
2064: PetscSectionGetFieldName(s, f, &name);
2065: PetscSectionSetFieldName(sNew, f, name);
2066: PetscSectionGetFieldComponents(s, f, &numComp);
2067: PetscSectionSetFieldComponents(sNew, f, numComp);
2068: }
2069: ISGetLocalSize(permutation, &numPoints);
2070: ISGetIndices(permutation, &perm);
2071: PetscSectionGetChart(s, &pStart, &pEnd);
2072: PetscSectionSetChart(sNew, pStart, pEnd);
2073: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
2074: for (p = pStart; p < pEnd; ++p) {
2075: PetscInt dof, cdof;
2077: PetscSectionGetDof(s, p, &dof);
2078: PetscSectionSetDof(sNew, perm[p], dof);
2079: PetscSectionGetConstraintDof(s, p, &cdof);
2080: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2081: for (f = 0; f < numFields; ++f) {
2082: PetscSectionGetFieldDof(s, p, f, &dof);
2083: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2084: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2085: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2086: }
2087: }
2088: PetscSectionSetUp(sNew);
2089: for (p = pStart; p < pEnd; ++p) {
2090: const PetscInt *cind;
2091: PetscInt cdof;
2093: PetscSectionGetConstraintDof(s, p, &cdof);
2094: if (cdof) {
2095: PetscSectionGetConstraintIndices(s, p, &cind);
2096: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2097: }
2098: for (f = 0; f < numFields; ++f) {
2099: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2100: if (cdof) {
2101: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2102: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2103: }
2104: }
2105: }
2106: ISRestoreIndices(permutation, &perm);
2107: *sectionNew = sNew;
2108: return(0);
2109: }
2111: /*@C
2112: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
2114: Collective
2116: Input Parameters:
2117: + sf - The SF
2118: - rootSection - Section defined on root space
2120: Output Parameters:
2121: + remoteOffsets - root offsets in leaf storage, or NULL
2122: - leafSection - Section defined on the leaf space
2124: Level: intermediate
2126: .seealso: PetscSFCreate()
2127: @*/
2128: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2129: {
2130: PetscSF embedSF;
2131: const PetscInt *ilocal, *indices;
2132: IS selected;
2133: PetscInt numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;
2137: PetscSectionGetNumFields(rootSection, &numFields);
2138: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2139: for (f = 0; f < numFields; ++f) {
2140: PetscInt numComp = 0;
2141: PetscSectionGetFieldComponents(rootSection, f, &numComp);
2142: PetscSectionSetFieldComponents(leafSection, f, numComp);
2143: }
2144: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2145: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2146: rpEnd = PetscMin(rpEnd,nroots);
2147: rpEnd = PetscMax(rpStart,rpEnd);
2148: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2149: ISGetIndices(selected, &indices);
2150: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2151: ISRestoreIndices(selected, &indices);
2152: ISDestroy(&selected);
2153: PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
2154: if (nleaves && ilocal) {
2155: for (i = 0; i < nleaves; ++i) {
2156: lpStart = PetscMin(lpStart, ilocal[i]);
2157: lpEnd = PetscMax(lpEnd, ilocal[i]);
2158: }
2159: ++lpEnd;
2160: } else {
2161: lpStart = 0;
2162: lpEnd = nleaves;
2163: }
2164: PetscSectionSetChart(leafSection, lpStart, lpEnd);
2165: /* Could fuse these at the cost of a copy and extra allocation */
2166: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2167: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2168: for (f = 0; f < numFields; ++f) {
2169: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2170: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2171: }
2172: if (remoteOffsets) {
2173: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2174: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2175: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2176: }
2177: PetscSFDestroy(&embedSF);
2178: PetscSectionSetUp(leafSection);
2179: return(0);
2180: }
2182: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2183: {
2184: PetscSF embedSF;
2185: const PetscInt *indices;
2186: IS selected;
2187: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2188: PetscErrorCode ierr;
2191: *remoteOffsets = NULL;
2192: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2193: if (numRoots < 0) return(0);
2194: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2195: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2196: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2197: ISGetIndices(selected, &indices);
2198: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2199: ISRestoreIndices(selected, &indices);
2200: ISDestroy(&selected);
2201: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2202: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2203: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2204: PetscSFDestroy(&embedSF);
2205: return(0);
2206: }
2208: /*@C
2209: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2211: Input Parameters:
2212: + sf - The SF
2213: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2214: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2215: - leafSection - Data layout of local points for incoming data (this is the distributed section)
2217: Output Parameters:
2218: - sectionSF - The new SF
2220: Note: Either rootSection or remoteOffsets can be specified
2222: Level: intermediate
2224: .seealso: PetscSFCreate()
2225: @*/
2226: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2227: {
2228: MPI_Comm comm;
2229: const PetscInt *localPoints;
2230: const PetscSFNode *remotePoints;
2231: PetscInt lpStart, lpEnd;
2232: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2233: PetscInt *localIndices;
2234: PetscSFNode *remoteIndices;
2235: PetscInt i, ind;
2236: PetscErrorCode ierr;
2244: PetscObjectGetComm((PetscObject)sf,&comm);
2245: PetscSFCreate(comm, sectionSF);
2246: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2247: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2248: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2249: if (numRoots < 0) return(0);
2250: for (i = 0; i < numPoints; ++i) {
2251: PetscInt localPoint = localPoints ? localPoints[i] : i;
2252: PetscInt dof;
2254: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2255: PetscSectionGetDof(leafSection, localPoint, &dof);
2256: numIndices += dof;
2257: }
2258: }
2259: PetscMalloc1(numIndices, &localIndices);
2260: PetscMalloc1(numIndices, &remoteIndices);
2261: /* Create new index graph */
2262: for (i = 0, ind = 0; i < numPoints; ++i) {
2263: PetscInt localPoint = localPoints ? localPoints[i] : i;
2264: PetscInt rank = remotePoints[i].rank;
2266: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2267: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2268: PetscInt loff, dof, d;
2270: PetscSectionGetOffset(leafSection, localPoint, &loff);
2271: PetscSectionGetDof(leafSection, localPoint, &dof);
2272: for (d = 0; d < dof; ++d, ++ind) {
2273: localIndices[ind] = loff+d;
2274: remoteIndices[ind].rank = rank;
2275: remoteIndices[ind].index = remoteOffset+d;
2276: }
2277: }
2278: }
2279: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2280: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2281: return(0);
2282: }
2284: /*@
2285: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2287: Input Parameters:
2288: + section - The PetscSection
2289: . obj - A PetscObject which serves as the key for this index
2290: . clSection - Section giving the size of the closure of each point
2291: - clPoints - IS giving the points in each closure
2293: Note: We compress out closure points with no dofs in this section
2295: Level: intermediate
2297: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2298: @*/
2299: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2300: {
2304: if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2305: section->clObj = obj;
2306: PetscSectionDestroy(§ion->clSection);
2307: ISDestroy(§ion->clPoints);
2308: section->clSection = clSection;
2309: section->clPoints = clPoints;
2310: return(0);
2311: }
2313: /*@
2314: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2316: Input Parameters:
2317: + section - The PetscSection
2318: - obj - A PetscObject which serves as the key for this index
2320: Output Parameters:
2321: + clSection - Section giving the size of the closure of each point
2322: - clPoints - IS giving the points in each closure
2324: Note: We compress out closure points with no dofs in this section
2326: Level: intermediate
2328: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2329: @*/
2330: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2331: {
2333: if (section->clObj == obj) {
2334: if (clSection) *clSection = section->clSection;
2335: if (clPoints) *clPoints = section->clPoints;
2336: } else {
2337: if (clSection) *clSection = NULL;
2338: if (clPoints) *clPoints = NULL;
2339: }
2340: return(0);
2341: }
2343: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2344: {
2345: PetscInt i;
2349: if (section->clObj != obj) {
2350: PetscSectionDestroy(§ion->clSection);
2351: ISDestroy(§ion->clPoints);
2352: }
2353: section->clObj = obj;
2354: PetscFree(section->clPerm);
2355: PetscFree(section->clInvPerm);
2356: section->clSize = clSize;
2357: if (mode == PETSC_COPY_VALUES) {
2358: PetscMalloc1(clSize, §ion->clPerm);
2359: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2360: PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2361: } else if (mode == PETSC_OWN_POINTER) {
2362: section->clPerm = clPerm;
2363: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2364: PetscMalloc1(clSize, §ion->clInvPerm);
2365: for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2366: return(0);
2367: }
2369: /*@
2370: PetscSectionSetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2372: Input Parameters:
2373: + section - The PetscSection
2374: . obj - A PetscObject which serves as the key for this index
2375: - perm - Permutation of the cell dof closure
2377: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2378: other points (like faces).
2380: Level: intermediate
2382: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2383: @*/
2384: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2385: {
2386: const PetscInt *clPerm = NULL;
2387: PetscInt clSize = 0;
2388: PetscErrorCode ierr;
2391: if (perm) {
2392: ISGetLocalSize(perm, &clSize);
2393: ISGetIndices(perm, &clPerm);
2394: }
2395: PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2396: if (perm) {ISRestoreIndices(perm, &clPerm);}
2397: return(0);
2398: }
2400: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2401: {
2403: if (section->clObj == obj) {
2404: if (size) *size = section->clSize;
2405: if (perm) *perm = section->clPerm;
2406: } else {
2407: if (size) *size = 0;
2408: if (perm) *perm = NULL;
2409: }
2410: return(0);
2411: }
2413: /*@
2414: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2416: Input Parameters:
2417: + section - The PetscSection
2418: - obj - A PetscObject which serves as the key for this index
2420: Output Parameter:
2421: . perm - The dof closure permutation
2423: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2424: other points (like faces).
2426: The user must destroy the IS that is returned.
2428: Level: intermediate
2430: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2431: @*/
2432: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2433: {
2434: const PetscInt *clPerm;
2435: PetscInt clSize;
2436: PetscErrorCode ierr;
2439: PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2440: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2441: return(0);
2442: }
2444: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2445: {
2447: if (section->clObj == obj) {
2448: if (size) *size = section->clSize;
2449: if (perm) *perm = section->clInvPerm;
2450: } else {
2451: if (size) *size = 0;
2452: if (perm) *perm = NULL;
2453: }
2454: return(0);
2455: }
2457: /*@
2458: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2460: Input Parameters:
2461: + section - The PetscSection
2462: - obj - A PetscObject which serves as the key for this index
2464: Output Parameters:
2465: + size - The dof closure size
2466: - perm - The dof closure permutation
2468: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2469: other points (like faces).
2471: The user must destroy the IS that is returned.
2473: Level: intermediate
2475: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2476: @*/
2477: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2478: {
2479: const PetscInt *clPerm;
2480: PetscInt clSize;
2481: PetscErrorCode ierr;
2484: PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2485: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2486: return(0);
2487: }
2489: /*@
2490: PetscSectionGetField - Get the subsection associated with a single field
2492: Input Parameters:
2493: + s - The PetscSection
2494: - field - The field number
2496: Output Parameter:
2497: . subs - The subsection for the given field
2499: Level: intermediate
2501: .seealso: PetscSectionSetNumFields()
2502: @*/
2503: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2504: {
2510: 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);
2511: PetscObjectReference((PetscObject) s->field[field]);
2512: *subs = s->field[field];
2513: return(0);
2514: }
2516: PetscClassId PETSC_SECTION_SYM_CLASSID;
2517: PetscFunctionList PetscSectionSymList = NULL;
2519: /*@
2520: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
2522: Collective on MPI_Comm
2524: Input Parameter:
2525: . comm - the MPI communicator
2527: Output Parameter:
2528: . sym - pointer to the new set of symmetries
2530: Level: developer
2532: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2533: @*/
2534: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2535: {
2540: ISInitializePackage();
2541: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2542: return(0);
2543: }
2545: /*@C
2546: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
2548: Collective on PetscSectionSym
2550: Input Parameters:
2551: + sym - The section symmetry object
2552: - method - The name of the section symmetry type
2554: Level: developer
2556: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2557: @*/
2558: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2559: {
2560: PetscErrorCode (*r)(PetscSectionSym);
2561: PetscBool match;
2566: PetscObjectTypeCompare((PetscObject) sym, method, &match);
2567: if (match) return(0);
2569: PetscFunctionListFind(PetscSectionSymList,method,&r);
2570: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2571: if (sym->ops->destroy) {
2572: (*sym->ops->destroy)(sym);
2573: sym->ops->destroy = NULL;
2574: }
2575: (*r)(sym);
2576: PetscObjectChangeTypeName((PetscObject)sym,method);
2577: return(0);
2578: }
2581: /*@C
2582: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
2584: Not Collective
2586: Input Parameter:
2587: . sym - The section symmetry
2589: Output Parameter:
2590: . type - The index set type name
2592: Level: developer
2594: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2595: @*/
2596: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2597: {
2601: *type = ((PetscObject)sym)->type_name;
2602: return(0);
2603: }
2605: /*@C
2606: PetscSectionSymRegister - Adds a new section symmetry implementation
2608: Not Collective
2610: Input Parameters:
2611: + name - The name of a new user-defined creation routine
2612: - create_func - The creation routine itself
2614: Notes:
2615: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
2617: Level: developer
2619: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2620: @*/
2621: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2622: {
2626: ISInitializePackage();
2627: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2628: return(0);
2629: }
2631: /*@
2632: PetscSectionSymDestroy - Destroys a section symmetry.
2634: Collective on PetscSectionSym
2636: Input Parameters:
2637: . sym - the section symmetry
2639: Level: developer
2641: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2642: @*/
2643: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2644: {
2645: SymWorkLink link,next;
2649: if (!*sym) return(0);
2651: if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2652: if ((*sym)->ops->destroy) {
2653: (*(*sym)->ops->destroy)(*sym);
2654: }
2655: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2656: for (link=(*sym)->workin; link; link=next) {
2657: next = link->next;
2658: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2659: PetscFree(link);
2660: }
2661: (*sym)->workin = NULL;
2662: PetscHeaderDestroy(sym);
2663: return(0);
2664: }
2666: /*@C
2667: PetscSectionSymView - Displays a section symmetry
2669: Collective on PetscSectionSym
2671: Input Parameters:
2672: + sym - the index set
2673: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
2675: Level: developer
2677: .seealso: PetscViewerASCIIOpen()
2678: @*/
2679: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2680: {
2685: if (!viewer) {
2686: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2687: }
2690: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2691: if (sym->ops->view) {
2692: (*sym->ops->view)(sym,viewer);
2693: }
2694: return(0);
2695: }
2697: /*@
2698: PetscSectionSetSym - Set the symmetries for the data referred to by the section
2700: Collective on PetscSection
2702: Input Parameters:
2703: + section - the section describing data layout
2704: - sym - the symmetry describing the affect of orientation on the access of the data
2706: Level: developer
2708: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2709: @*/
2710: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2711: {
2716: PetscSectionSymDestroy(&(section->sym));
2717: if (sym) {
2720: PetscObjectReference((PetscObject) sym);
2721: }
2722: section->sym = sym;
2723: return(0);
2724: }
2726: /*@
2727: PetscSectionGetSym - Get the symmetries for the data referred to by the section
2729: Not collective
2731: Input Parameters:
2732: . section - the section describing data layout
2734: Output Parameters:
2735: . sym - the symmetry describing the affect of orientation on the access of the data
2737: Level: developer
2739: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2740: @*/
2741: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2742: {
2745: *sym = section->sym;
2746: return(0);
2747: }
2749: /*@
2750: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
2752: Collective on PetscSection
2754: Input Parameters:
2755: + section - the section describing data layout
2756: . field - the field number
2757: - sym - the symmetry describing the affect of orientation on the access of the data
2759: Level: developer
2761: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2762: @*/
2763: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2764: {
2769: 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);
2770: PetscSectionSetSym(section->field[field],sym);
2771: return(0);
2772: }
2774: /*@
2775: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
2777: Collective on PetscSection
2779: Input Parameters:
2780: + section - the section describing data layout
2781: - field - the field number
2783: Output Parameters:
2784: . sym - the symmetry describing the affect of orientation on the access of the data
2786: Level: developer
2788: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2789: @*/
2790: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2791: {
2794: 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);
2795: *sym = section->field[field]->sym;
2796: return(0);
2797: }
2799: /*@C
2800: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
2802: Not collective
2804: Input Parameters:
2805: + section - the section
2806: . numPoints - the number of points
2807: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2808: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2809: context, see DMPlexGetConeOrientation()).
2811: Output Parameter:
2812: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2813: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2814: identity).
2816: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2817: .vb
2818: const PetscInt **perms;
2819: const PetscScalar **rots;
2820: PetscInt lOffset;
2822: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2823: for (i = 0, lOffset = 0; i < numPoints; i++) {
2824: PetscInt point = points[2*i], dof, sOffset;
2825: const PetscInt *perm = perms ? perms[i] : NULL;
2826: const PetscScalar *rot = rots ? rots[i] : NULL;
2828: PetscSectionGetDof(section,point,&dof);
2829: PetscSectionGetOffset(section,point,&sOffset);
2831: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
2832: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
2833: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
2834: lOffset += dof;
2835: }
2836: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2837: .ve
2839: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2840: .vb
2841: const PetscInt **perms;
2842: const PetscScalar **rots;
2843: PetscInt lOffset;
2845: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2846: for (i = 0, lOffset = 0; i < numPoints; i++) {
2847: PetscInt point = points[2*i], dof, sOffset;
2848: const PetscInt *perm = perms ? perms[i] : NULL;
2849: const PetscScalar *rot = rots ? rots[i] : NULL;
2851: PetscSectionGetDof(section,point,&dof);
2852: PetscSectionGetOffset(section,point,&sOff);
2854: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2855: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
2856: offset += dof;
2857: }
2858: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2859: .ve
2861: Level: developer
2863: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2864: @*/
2865: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2866: {
2867: PetscSectionSym sym;
2868: PetscErrorCode ierr;
2873: if (perms) *perms = NULL;
2874: if (rots) *rots = NULL;
2875: sym = section->sym;
2876: if (sym && (perms || rots)) {
2877: SymWorkLink link;
2879: if (sym->workin) {
2880: link = sym->workin;
2881: sym->workin = sym->workin->next;
2882: } else {
2883: PetscNewLog(sym,&link);
2884: }
2885: if (numPoints > link->numPoints) {
2886: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2887: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
2888: link->numPoints = numPoints;
2889: }
2890: link->next = sym->workout;
2891: sym->workout = link;
2892: PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2893: PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2894: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2895: if (perms) *perms = link->perms;
2896: if (rots) *rots = link->rots;
2897: }
2898: return(0);
2899: }
2901: /*@C
2902: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
2904: Not collective
2906: Input Parameters:
2907: + section - the section
2908: . numPoints - the number of points
2909: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2910: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2911: context, see DMPlexGetConeOrientation()).
2913: Output Parameter:
2914: + perms - The permutations for the given orientations: set to NULL at conclusion
2915: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
2917: Level: developer
2919: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2920: @*/
2921: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2922: {
2923: PetscSectionSym sym;
2927: sym = section->sym;
2928: if (sym && (perms || rots)) {
2929: SymWorkLink *p,link;
2931: for (p=&sym->workout; (link=*p); p=&link->next) {
2932: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
2933: *p = link->next;
2934: link->next = sym->workin;
2935: sym->workin = link;
2936: if (perms) *perms = NULL;
2937: if (rots) *rots = NULL;
2938: return(0);
2939: }
2940: }
2941: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
2942: }
2943: return(0);
2944: }
2946: /*@C
2947: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
2949: Not collective
2951: Input Parameters:
2952: + section - the section
2953: . field - the field of the section
2954: . numPoints - the number of points
2955: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2956: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2957: context, see DMPlexGetConeOrientation()).
2959: Output Parameter:
2960: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2961: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2962: identity).
2964: Level: developer
2966: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
2967: @*/
2968: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2969: {
2974: 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);
2975: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
2976: return(0);
2977: }
2979: /*@C
2980: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
2982: Not collective
2984: Input Parameters:
2985: + section - the section
2986: . field - the field number
2987: . numPoints - the number of points
2988: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2989: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2990: context, see DMPlexGetConeOrientation()).
2992: Output Parameter:
2993: + perms - The permutations for the given orientations: set to NULL at conclusion
2994: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
2996: Level: developer
2998: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2999: @*/
3000: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3001: {
3006: 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);
3007: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3008: return(0);
3009: }
3011: /*@
3012: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3014: Not collective
3016: Input Parameter:
3017: . s - the global PetscSection
3019: Output Parameters:
3020: . flg - the flag
3022: Level: developer
3024: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3025: @*/
3026: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3027: {
3030: *flg = s->useFieldOff;
3031: return(0);
3032: }
3034: /*@
3035: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3037: Not collective
3039: Input Parameters:
3040: + s - the global PetscSection
3041: - flg - the flag
3043: Level: developer
3045: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3046: @*/
3047: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3048: {
3051: s->useFieldOff = flg;
3052: return(0);
3053: }