Actual source code: vsectionis.c
petsc-3.8.4 2018-03-24
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: Typical calling sequence
23: $ PetscSectionCreate(MPI_Comm,PetscSection *);
24: $ PetscSectionSetNumFields(PetscSection, numFields);
25: $ PetscSectionSetChart(PetscSection,low,high);
26: $ PetscSectionSetDof(PetscSection,point,numdof);
27: $ PetscSectionSetUp(PetscSection);
28: $ PetscSectionGetOffset(PetscSection,point,PetscInt *);
29: $ PetscSectionDestroy(PetscSection);
31: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
32: recommended they not be used in user codes unless you really gain something in their use.
34: .seealso: PetscSection, PetscSectionDestroy()
35: @*/
36: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
37: {
42: ISInitializePackage();
44: PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
46: (*s)->pStart = -1;
47: (*s)->pEnd = -1;
48: (*s)->perm = NULL;
49: (*s)->maxDof = 0;
50: (*s)->atlasDof = NULL;
51: (*s)->atlasOff = NULL;
52: (*s)->bc = NULL;
53: (*s)->bcIndices = NULL;
54: (*s)->setup = PETSC_FALSE;
55: (*s)->numFields = 0;
56: (*s)->fieldNames = NULL;
57: (*s)->field = NULL;
58: (*s)->clObj = NULL;
59: (*s)->clSection = NULL;
60: (*s)->clPoints = NULL;
61: (*s)->clSize = 0;
62: (*s)->clPerm = NULL;
63: (*s)->clInvPerm = NULL;
64: return(0);
65: }
67: /*@
68: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
70: Collective on MPI_Comm
72: Input Parameter:
73: . section - the PetscSection
75: Output Parameter:
76: . newSection - the copy
78: Level: developer
80: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
81: @*/
82: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
83: {
84: IS perm;
85: PetscInt numFields, f, pStart, pEnd, p;
89: PetscSectionGetNumFields(section, &numFields);
90: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
91: for (f = 0; f < numFields; ++f) {
92: const char *name = NULL;
93: PetscInt numComp = 0;
95: PetscSectionGetFieldName(section, f, &name);
96: PetscSectionSetFieldName(newSection, f, name);
97: PetscSectionGetFieldComponents(section, f, &numComp);
98: PetscSectionSetFieldComponents(newSection, f, numComp);
99: }
100: PetscSectionGetChart(section, &pStart, &pEnd);
101: PetscSectionSetChart(newSection, pStart, pEnd);
102: PetscSectionGetPermutation(section, &perm);
103: PetscSectionSetPermutation(newSection, perm);
104: for (p = pStart; p < pEnd; ++p) {
105: PetscInt dof, cdof, fcdof = 0;
107: PetscSectionGetDof(section, p, &dof);
108: PetscSectionSetDof(newSection, p, dof);
109: PetscSectionGetConstraintDof(section, p, &cdof);
110: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
111: for (f = 0; f < numFields; ++f) {
112: PetscSectionGetFieldDof(section, p, f, &dof);
113: PetscSectionSetFieldDof(newSection, p, f, dof);
114: if (cdof) {
115: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
116: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
117: }
118: }
119: }
120: PetscSectionSetUp(newSection);
121: for (p = pStart; p < pEnd; ++p) {
122: PetscInt off, cdof, fcdof = 0;
123: const PetscInt *cInd;
125: /* Must set offsets in case they do not agree with the prefix sums */
126: PetscSectionGetOffset(section, p, &off);
127: PetscSectionSetOffset(newSection, p, off);
128: PetscSectionGetConstraintDof(section, p, &cdof);
129: if (cdof) {
130: PetscSectionGetConstraintIndices(section, p, &cInd);
131: PetscSectionSetConstraintIndices(newSection, p, cInd);
132: for (f = 0; f < numFields; ++f) {
133: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
134: if (fcdof) {
135: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
136: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
137: }
138: }
139: }
140: }
141: return(0);
142: }
144: /*@
145: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
147: Collective on MPI_Comm
149: Input Parameter:
150: . section - the PetscSection
152: Output Parameter:
153: . newSection - the copy
155: Level: developer
157: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
158: @*/
159: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
160: {
164: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
165: PetscSectionCopy(section, *newSection);
166: return(0);
167: }
169: /*@
170: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
172: Not collective
174: Input Parameter:
175: . s - the PetscSection
177: Output Parameter:
178: . numFields - the number of fields defined, or 0 if none were defined
180: Level: intermediate
182: .seealso: PetscSectionSetNumFields()
183: @*/
184: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
185: {
188: *numFields = s->numFields;
189: return(0);
190: }
192: /*@
193: PetscSectionSetNumFields - Sets the number of fields.
195: Not collective
197: Input Parameters:
198: + s - the PetscSection
199: - numFields - the number of fields
201: Level: intermediate
203: .seealso: PetscSectionGetNumFields()
204: @*/
205: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
206: {
207: PetscInt f;
211: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %d must be positive", numFields);
212: PetscSectionReset(s);
214: s->numFields = numFields;
215: PetscMalloc1(s->numFields, &s->numFieldComponents);
216: PetscMalloc1(s->numFields, &s->fieldNames);
217: PetscMalloc1(s->numFields, &s->field);
218: for (f = 0; f < s->numFields; ++f) {
219: char name[64];
221: s->numFieldComponents[f] = 1;
223: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
224: PetscSNPrintf(name, 64, "Field_%D", f);
225: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
226: }
227: return(0);
228: }
230: /*@C
231: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
233: Not Collective
235: Input Parameters:
236: + s - the PetscSection
237: - field - the field number
239: Output Parameter:
240: . fieldName - the field name
242: Level: developer
244: .seealso: PetscSectionSetFieldName()
245: @*/
246: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
247: {
250: 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);
251: *fieldName = s->fieldNames[field];
252: return(0);
253: }
255: /*@C
256: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
258: Not Collective
260: Input Parameters:
261: + s - the PetscSection
262: . field - the field number
263: - fieldName - the field name
265: Level: developer
267: .seealso: PetscSectionGetFieldName()
268: @*/
269: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
270: {
275: 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);
276: PetscFree(s->fieldNames[field]);
277: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
278: return(0);
279: }
281: /*@
282: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
284: Not collective
286: Input Parameters:
287: + s - the PetscSection
288: - field - the field number
290: Output Parameter:
291: . numComp - the number of field components
293: Level: intermediate
295: .seealso: PetscSectionSetNumFieldComponents(), PetscSectionGetNumFields()
296: @*/
297: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
298: {
301: 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);
302: *numComp = s->numFieldComponents[field];
303: return(0);
304: }
306: /*@
307: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
309: Not collective
311: Input Parameters:
312: + s - the PetscSection
313: . field - the field number
314: - numComp - the number of field components
316: Level: intermediate
318: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
319: @*/
320: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
321: {
323: 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);
324: s->numFieldComponents[field] = numComp;
325: return(0);
326: }
328: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
329: {
333: if (!s->bc) {
334: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
335: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
336: }
337: return(0);
338: }
340: /*@
341: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
343: Not collective
345: Input Parameter:
346: . s - the PetscSection
348: Output Parameters:
349: + pStart - the first point
350: - pEnd - one past the last point
352: Level: intermediate
354: .seealso: PetscSectionSetChart(), PetscSectionCreate()
355: @*/
356: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
357: {
360: if (pStart) *pStart = s->pStart;
361: if (pEnd) *pEnd = s->pEnd;
362: return(0);
363: }
365: /*@
366: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
368: Not collective
370: Input Parameters:
371: + s - the PetscSection
372: . pStart - the first point
373: - pEnd - one past the last point
375: Level: intermediate
377: .seealso: PetscSectionGetChart(), PetscSectionCreate()
378: @*/
379: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
380: {
381: PetscInt f;
386: /* Cannot Reset() because it destroys field information */
387: s->setup = PETSC_FALSE;
388: PetscSectionDestroy(&s->bc);
389: PetscFree(s->bcIndices);
390: PetscFree2(s->atlasDof, s->atlasOff);
392: s->pStart = pStart;
393: s->pEnd = pEnd;
394: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
395: PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
396: for (f = 0; f < s->numFields; ++f) {
397: PetscSectionSetChart(s->field[f], pStart, pEnd);
398: }
399: return(0);
400: }
402: /*@
403: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
405: Not collective
407: Input Parameter:
408: . s - the PetscSection
410: Output Parameters:
411: . perm - The permutation as an IS
413: Level: intermediate
415: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
416: @*/
417: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
418: {
422: return(0);
423: }
425: /*@
426: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
428: Not collective
430: Input Parameters:
431: + s - the PetscSection
432: - perm - the permutation of points
434: Level: intermediate
436: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
437: @*/
438: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
439: {
443: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
444: if (s->perm != perm) {
445: ISDestroy(&s->perm);
446: s->perm = perm;
447: PetscObjectReference((PetscObject) s->perm);
448: }
449: return(0);
450: }
452: /*@
453: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
455: Not collective
457: Input Parameters:
458: + s - the PetscSection
459: - point - the point
461: Output Parameter:
462: . numDof - the number of dof
464: Level: intermediate
466: .seealso: PetscSectionSetDof(), PetscSectionCreate()
467: @*/
468: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
469: {
471: #if defined(PETSC_USE_DEBUG)
472: 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);
473: #endif
474: *numDof = s->atlasDof[point - s->pStart];
475: return(0);
476: }
478: /*@
479: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
481: Not collective
483: Input Parameters:
484: + s - the PetscSection
485: . point - the point
486: - numDof - the number of dof
488: Level: intermediate
490: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
491: @*/
492: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
493: {
495: 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);
496: s->atlasDof[point - s->pStart] = numDof;
497: return(0);
498: }
500: /*@
501: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
503: Not collective
505: Input Parameters:
506: + s - the PetscSection
507: . point - the point
508: - numDof - the number of additional dof
510: Level: intermediate
512: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
513: @*/
514: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
515: {
517: 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);
518: s->atlasDof[point - s->pStart] += numDof;
519: return(0);
520: }
522: /*@
523: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
525: Not collective
527: Input Parameters:
528: + s - the PetscSection
529: . point - the point
530: - field - the field
532: Output Parameter:
533: . numDof - the number of dof
535: Level: intermediate
537: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
538: @*/
539: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
540: {
544: 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);
545: PetscSectionGetDof(s->field[field], point, numDof);
546: return(0);
547: }
549: /*@
550: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
552: Not collective
554: Input Parameters:
555: + s - the PetscSection
556: . point - the point
557: . field - the field
558: - numDof - the number of dof
560: Level: intermediate
562: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
563: @*/
564: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
565: {
569: 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);
570: PetscSectionSetDof(s->field[field], point, numDof);
571: return(0);
572: }
574: /*@
575: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
577: Not collective
579: Input Parameters:
580: + s - the PetscSection
581: . point - the point
582: . field - the field
583: - numDof - the number of dof
585: Level: intermediate
587: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
588: @*/
589: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
590: {
594: 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);
595: PetscSectionAddDof(s->field[field], point, numDof);
596: return(0);
597: }
599: /*@
600: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
602: Not collective
604: Input Parameters:
605: + s - the PetscSection
606: - point - the point
608: Output Parameter:
609: . numDof - the number of dof which are fixed by constraints
611: Level: intermediate
613: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
614: @*/
615: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
616: {
620: if (s->bc) {
621: PetscSectionGetDof(s->bc, point, numDof);
622: } else *numDof = 0;
623: return(0);
624: }
626: /*@
627: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
629: Not collective
631: Input Parameters:
632: + s - the PetscSection
633: . point - the point
634: - numDof - the number of dof which are fixed by constraints
636: Level: intermediate
638: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
639: @*/
640: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
641: {
645: if (numDof) {
646: PetscSectionCheckConstraints_Static(s);
647: PetscSectionSetDof(s->bc, point, numDof);
648: }
649: return(0);
650: }
652: /*@
653: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
655: Not collective
657: Input Parameters:
658: + s - the PetscSection
659: . point - the point
660: - numDof - the number of additional dof which are fixed by constraints
662: Level: intermediate
664: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
665: @*/
666: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
667: {
671: if (numDof) {
672: PetscSectionCheckConstraints_Static(s);
673: PetscSectionAddDof(s->bc, point, numDof);
674: }
675: return(0);
676: }
678: /*@
679: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
681: Not collective
683: Input Parameters:
684: + s - the PetscSection
685: . point - the point
686: - field - the field
688: Output Parameter:
689: . numDof - the number of dof which are fixed by constraints
691: Level: intermediate
693: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
694: @*/
695: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
696: {
700: 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);
701: PetscSectionGetConstraintDof(s->field[field], point, numDof);
702: return(0);
703: }
705: /*@
706: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
708: Not collective
710: Input Parameters:
711: + s - the PetscSection
712: . point - the point
713: . field - the field
714: - numDof - the number of dof which are fixed by constraints
716: Level: intermediate
718: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
719: @*/
720: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
721: {
725: 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);
726: PetscSectionSetConstraintDof(s->field[field], point, numDof);
727: return(0);
728: }
730: /*@
731: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
733: Not collective
735: Input Parameters:
736: + s - the PetscSection
737: . point - the point
738: . field - the field
739: - numDof - the number of additional dof which are fixed by constraints
741: Level: intermediate
743: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
744: @*/
745: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
746: {
750: 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);
751: PetscSectionAddConstraintDof(s->field[field], point, numDof);
752: return(0);
753: }
755: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
756: {
760: if (s->bc) {
761: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
763: PetscSectionSetUp(s->bc);
764: PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
765: }
766: return(0);
767: }
769: /*@
770: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
772: Not collective
774: Input Parameter:
775: . s - the PetscSection
777: Level: intermediate
779: .seealso: PetscSectionCreate()
780: @*/
781: PetscErrorCode PetscSectionSetUp(PetscSection s)
782: {
783: const PetscInt *pind = NULL;
784: PetscInt offset = 0, p, f;
785: PetscErrorCode ierr;
788: if (s->setup) return(0);
789: s->setup = PETSC_TRUE;
790: if (s->perm) {ISGetIndices(s->perm, &pind);}
791: for (p = 0; p < s->pEnd - s->pStart; ++p) {
792: const PetscInt q = pind ? pind[p] : p;
794: s->atlasOff[q] = offset;
795: offset += s->atlasDof[q];
796: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
797: }
798: PetscSectionSetUpBC(s);
799: /* Assume that all fields have the same chart */
800: for (p = 0; p < s->pEnd - s->pStart; ++p) {
801: const PetscInt q = pind ? pind[p] : p;
802: PetscInt off = s->atlasOff[q];
804: for (f = 0; f < s->numFields; ++f) {
805: PetscSection sf = s->field[f];
807: sf->atlasOff[q] = off;
808: off += sf->atlasDof[q];
809: }
810: }
811: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
812: for (f = 0; f < s->numFields; ++f) {
813: PetscSectionSetUpBC(s->field[f]);
814: }
815: return(0);
816: }
818: /*@
819: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
821: Not collective
823: Input Parameters:
824: . s - the PetscSection
826: Output Parameter:
827: . maxDof - the maximum dof
829: Level: intermediate
831: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
832: @*/
833: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
834: {
836: *maxDof = s->maxDof;
837: return(0);
838: }
840: /*@
841: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
843: Not collective
845: Input Parameters:
846: + s - the PetscSection
847: - size - the allocated size
849: Output Parameter:
850: . size - the size of an array which can hold all the dofs
852: Level: intermediate
854: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
855: @*/
856: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
857: {
858: PetscInt p, n = 0;
861: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
862: *size = n;
863: return(0);
864: }
866: /*@
867: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
869: Not collective
871: Input Parameters:
872: + s - the PetscSection
873: - point - the point
875: Output Parameter:
876: . size - the size of an array which can hold all unconstrained dofs
878: Level: intermediate
880: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
881: @*/
882: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
883: {
884: PetscInt p, n = 0;
887: for (p = 0; p < s->pEnd - s->pStart; ++p) {
888: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
889: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
890: }
891: *size = n;
892: return(0);
893: }
895: /*@
896: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
897: the local section and an SF describing the section point overlap.
899: Input Parameters:
900: + s - The PetscSection for the local field layout
901: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
902: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
903: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
905: Output Parameter:
906: . gsection - The PetscSection for the global field layout
908: Note: This gives negative sizes and offsets to points not owned by this process
910: Level: developer
912: .seealso: PetscSectionCreate()
913: @*/
914: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
915: {
916: const PetscInt *pind = NULL;
917: PetscInt *recv = NULL, *neg = NULL;
918: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
919: PetscErrorCode ierr;
922: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
923: PetscSectionGetChart(s, &pStart, &pEnd);
924: PetscSectionSetChart(*gsection, pStart, pEnd);
925: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
926: nlocal = nroots; /* The local/leaf space matches global/root space */
927: /* Must allocate for all points visible to SF, which may be more than this section */
928: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
929: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %d < pEnd %d", nroots, pEnd);
930: PetscMalloc2(nroots,&neg,nlocal,&recv);
931: PetscMemzero(neg,nroots*sizeof(PetscInt));
932: }
933: /* Mark all local points with negative dof */
934: for (p = pStart; p < pEnd; ++p) {
935: PetscSectionGetDof(s, p, &dof);
936: PetscSectionSetDof(*gsection, p, dof);
937: PetscSectionGetConstraintDof(s, p, &cdof);
938: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
939: if (neg) neg[p] = -(dof+1);
940: }
941: PetscSectionSetUpBC(*gsection);
942: if (nroots >= 0) {
943: PetscMemzero(recv,nlocal*sizeof(PetscInt));
944: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
945: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
946: for (p = pStart; p < pEnd; ++p) {
947: if (recv[p] < 0) {
948: (*gsection)->atlasDof[p-pStart] = recv[p];
949: PetscSectionGetDof(s, p, &dof);
950: 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);
951: }
952: }
953: }
954: /* Calculate new sizes, get proccess offset, and calculate point offsets */
955: if (s->perm) {ISGetIndices(s->perm, &pind);}
956: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
957: const PetscInt q = pind ? pind[p] : p;
959: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
960: (*gsection)->atlasOff[q] = off;
961: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
962: }
963: if (!localOffsets) {
964: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
965: globalOff -= off;
966: }
967: for (p = pStart, off = 0; p < pEnd; ++p) {
968: (*gsection)->atlasOff[p-pStart] += globalOff;
969: if (neg) neg[p] = -((*gsection)->atlasOff[p-pStart]+1);
970: }
971: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
972: /* Put in negative offsets for ghost points */
973: if (nroots >= 0) {
974: PetscMemzero(recv,nlocal*sizeof(PetscInt));
975: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
976: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
977: for (p = pStart; p < pEnd; ++p) {
978: if (recv[p] < 0) (*gsection)->atlasOff[p-pStart] = recv[p];
979: }
980: }
981: PetscFree2(neg,recv);
982: PetscSectionViewFromOptions(*gsection,NULL,"-global_section_view");
983: return(0);
984: }
986: /*@
987: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
988: the local section and an SF describing the section point overlap.
990: Input Parameters:
991: + s - The PetscSection for the local field layout
992: . sf - The SF describing parallel layout of the section points
993: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
994: . numExcludes - The number of exclusion ranges
995: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
997: Output Parameter:
998: . gsection - The PetscSection for the global field layout
1000: Note: This gives negative sizes and offsets to points not owned by this process
1002: Level: developer
1004: .seealso: PetscSectionCreate()
1005: @*/
1006: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1007: {
1008: const PetscInt *pind = NULL;
1009: PetscInt *neg = NULL, *tmpOff = NULL;
1010: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1011: PetscErrorCode ierr;
1014: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1015: PetscSectionGetChart(s, &pStart, &pEnd);
1016: PetscSectionSetChart(*gsection, pStart, pEnd);
1017: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1018: if (nroots >= 0) {
1019: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %d < %d section size", nroots, pEnd-pStart);
1020: PetscCalloc1(nroots, &neg);
1021: if (nroots > pEnd-pStart) {
1022: PetscCalloc1(nroots, &tmpOff);
1023: } else {
1024: tmpOff = &(*gsection)->atlasDof[-pStart];
1025: }
1026: }
1027: /* Mark ghost points with negative dof */
1028: for (p = pStart; p < pEnd; ++p) {
1029: for (e = 0; e < numExcludes; ++e) {
1030: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1031: PetscSectionSetDof(*gsection, p, 0);
1032: break;
1033: }
1034: }
1035: if (e < numExcludes) continue;
1036: PetscSectionGetDof(s, p, &dof);
1037: PetscSectionSetDof(*gsection, p, dof);
1038: PetscSectionGetConstraintDof(s, p, &cdof);
1039: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1040: if (neg) neg[p] = -(dof+1);
1041: }
1042: PetscSectionSetUpBC(*gsection);
1043: if (nroots >= 0) {
1044: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1045: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1046: if (nroots > pEnd - pStart) {
1047: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1048: }
1049: }
1050: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1051: if (s->perm) {ISGetIndices(s->perm, &pind);}
1052: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1053: const PetscInt q = pind ? pind[p] : p;
1055: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1056: (*gsection)->atlasOff[q] = off;
1057: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1058: }
1059: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1060: globalOff -= off;
1061: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1062: (*gsection)->atlasOff[p] += globalOff;
1063: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1064: }
1065: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1066: /* Put in negative offsets for ghost points */
1067: if (nroots >= 0) {
1068: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1069: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1070: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1071: if (nroots > pEnd - pStart) {
1072: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1073: }
1074: }
1075: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1076: PetscFree(neg);
1077: return(0);
1078: }
1080: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1081: {
1082: PetscInt pStart, pEnd, p, localSize = 0;
1086: PetscSectionGetChart(s, &pStart, &pEnd);
1087: for (p = pStart; p < pEnd; ++p) {
1088: PetscInt dof;
1090: PetscSectionGetDof(s, p, &dof);
1091: if (dof > 0) ++localSize;
1092: }
1093: PetscLayoutCreate(comm, layout);
1094: PetscLayoutSetLocalSize(*layout, localSize);
1095: PetscLayoutSetBlockSize(*layout, 1);
1096: PetscLayoutSetUp(*layout);
1097: return(0);
1098: }
1100: /*@
1101: PetscSectionGetValueLayout - Get the PetscLayout associated with a section, usually the default global section.
1103: Input Parameters:
1104: + comm - The MPI_Comm
1105: - s - The PetscSection
1107: Output Parameter:
1108: . layout - The layout for the section
1110: Level: developer
1112: .seealso: PetscSectionCreate()
1113: @*/
1114: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1115: {
1116: PetscInt pStart, pEnd, p, localSize = 0;
1120: PetscSectionGetChart(s, &pStart, &pEnd);
1121: for (p = pStart; p < pEnd; ++p) {
1122: PetscInt dof,cdof;
1124: PetscSectionGetDof(s, p, &dof);
1125: PetscSectionGetConstraintDof(s, p, &cdof);
1126: if (dof-cdof > 0) localSize += dof-cdof;
1127: }
1128: PetscLayoutCreate(comm, layout);
1129: PetscLayoutSetLocalSize(*layout, localSize);
1130: PetscLayoutSetBlockSize(*layout, 1);
1131: PetscLayoutSetUp(*layout);
1132: return(0);
1133: }
1135: /*@
1136: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1138: Not collective
1140: Input Parameters:
1141: + s - the PetscSection
1142: - point - the point
1144: Output Parameter:
1145: . offset - the offset
1147: Level: intermediate
1149: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1150: @*/
1151: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1152: {
1154: #if defined(PETSC_USE_DEBUG)
1155: 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);
1156: #endif
1157: *offset = s->atlasOff[point - s->pStart];
1158: return(0);
1159: }
1161: /*@
1162: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1164: Not collective
1166: Input Parameters:
1167: + s - the PetscSection
1168: . point - the point
1169: - offset - the offset
1171: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1173: Level: intermediate
1175: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1176: @*/
1177: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1178: {
1180: 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);
1181: s->atlasOff[point - s->pStart] = offset;
1182: return(0);
1183: }
1185: /*@
1186: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1188: Not collective
1190: Input Parameters:
1191: + s - the PetscSection
1192: . point - the point
1193: - field - the field
1195: Output Parameter:
1196: . offset - the offset
1198: Level: intermediate
1200: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1201: @*/
1202: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1203: {
1207: 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);
1208: PetscSectionGetOffset(s->field[field], point, offset);
1209: return(0);
1210: }
1212: /*@
1213: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1215: Not collective
1217: Input Parameters:
1218: + s - the PetscSection
1219: . point - the point
1220: . field - the field
1221: - offset - the offset
1223: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1225: Level: intermediate
1227: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1228: @*/
1229: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1230: {
1234: 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);
1235: PetscSectionSetOffset(s->field[field], point, offset);
1236: return(0);
1237: }
1239: /* This gives the offset on a point of the field, ignoring constraints */
1240: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1241: {
1242: PetscInt off, foff;
1246: 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);
1247: PetscSectionGetOffset(s, point, &off);
1248: PetscSectionGetOffset(s->field[field], point, &foff);
1249: *offset = foff - off;
1250: return(0);
1251: }
1253: /*@
1254: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1256: Not collective
1258: Input Parameter:
1259: . s - the PetscSection
1261: Output Parameters:
1262: + start - the minimum offset
1263: - end - one more than the maximum offset
1265: Level: intermediate
1267: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1268: @*/
1269: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1270: {
1271: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1275: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1276: PetscSectionGetChart(s, &pStart, &pEnd);
1277: for (p = 0; p < pEnd-pStart; ++p) {
1278: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1280: if (off >= 0) {
1281: os = PetscMin(os, off);
1282: oe = PetscMax(oe, off+dof);
1283: }
1284: }
1285: if (start) *start = os;
1286: if (end) *end = oe;
1287: return(0);
1288: }
1290: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, PetscInt fields[], PetscSection *subs)
1291: {
1292: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1296: if (!numFields) return(0);
1297: PetscSectionGetNumFields(s, &nF);
1298: if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of fields %d", numFields, nF);
1299: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1300: PetscSectionSetNumFields(*subs, numFields);
1301: for (f = 0; f < numFields; ++f) {
1302: const char *name = NULL;
1303: PetscInt numComp = 0;
1305: PetscSectionGetFieldName(s, fields[f], &name);
1306: PetscSectionSetFieldName(*subs, f, name);
1307: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1308: PetscSectionSetFieldComponents(*subs, f, numComp);
1309: }
1310: PetscSectionGetChart(s, &pStart, &pEnd);
1311: PetscSectionSetChart(*subs, pStart, pEnd);
1312: for (p = pStart; p < pEnd; ++p) {
1313: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1315: for (f = 0; f < numFields; ++f) {
1316: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1317: PetscSectionSetFieldDof(*subs, p, f, fdof);
1318: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1319: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1320: dof += fdof;
1321: cdof += cfdof;
1322: }
1323: PetscSectionSetDof(*subs, p, dof);
1324: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1325: maxCdof = PetscMax(cdof, maxCdof);
1326: }
1327: PetscSectionSetUp(*subs);
1328: if (maxCdof) {
1329: PetscInt *indices;
1331: PetscMalloc1(maxCdof, &indices);
1332: for (p = pStart; p < pEnd; ++p) {
1333: PetscInt cdof;
1335: PetscSectionGetConstraintDof(*subs, p, &cdof);
1336: if (cdof) {
1337: const PetscInt *oldIndices = NULL;
1338: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1340: for (f = 0; f < numFields; ++f) {
1341: PetscInt oldFoff = 0, oldf;
1343: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1344: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1345: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1346: /* This can be sped up if we assume sorted fields */
1347: for (oldf = 0; oldf < fields[f]; ++oldf) {
1348: PetscInt oldfdof = 0;
1349: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1350: oldFoff += oldfdof;
1351: }
1352: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1353: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1354: numConst += cfdof;
1355: fOff += fdof;
1356: }
1357: PetscSectionSetConstraintIndices(*subs, p, indices);
1358: }
1359: }
1360: PetscFree(indices);
1361: }
1362: return(0);
1363: }
1365: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1366: {
1367: const PetscInt *points = NULL, *indices = NULL;
1368: PetscInt numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;
1372: PetscSectionGetNumFields(s, &numFields);
1373: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1374: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1375: for (f = 0; f < numFields; ++f) {
1376: const char *name = NULL;
1377: PetscInt numComp = 0;
1379: PetscSectionGetFieldName(s, f, &name);
1380: PetscSectionSetFieldName(*subs, f, name);
1381: PetscSectionGetFieldComponents(s, f, &numComp);
1382: PetscSectionSetFieldComponents(*subs, f, numComp);
1383: }
1384: /* For right now, we do not try to squeeze the subchart */
1385: if (subpointMap) {
1386: ISGetSize(subpointMap, &numSubpoints);
1387: ISGetIndices(subpointMap, &points);
1388: }
1389: PetscSectionGetChart(s, &pStart, &pEnd);
1390: PetscSectionSetChart(*subs, 0, numSubpoints);
1391: for (p = pStart; p < pEnd; ++p) {
1392: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1394: PetscFindInt(p, numSubpoints, points, &subp);
1395: if (subp < 0) continue;
1396: for (f = 0; f < numFields; ++f) {
1397: PetscSectionGetFieldDof(s, p, f, &fdof);
1398: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1399: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1400: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1401: }
1402: PetscSectionGetDof(s, p, &dof);
1403: PetscSectionSetDof(*subs, subp, dof);
1404: PetscSectionGetConstraintDof(s, p, &cdof);
1405: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1406: }
1407: PetscSectionSetUp(*subs);
1408: /* Change offsets to original offsets */
1409: for (p = pStart; p < pEnd; ++p) {
1410: PetscInt off, foff = 0;
1412: PetscFindInt(p, numSubpoints, points, &subp);
1413: if (subp < 0) continue;
1414: for (f = 0; f < numFields; ++f) {
1415: PetscSectionGetFieldOffset(s, p, f, &foff);
1416: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1417: }
1418: PetscSectionGetOffset(s, p, &off);
1419: PetscSectionSetOffset(*subs, subp, off);
1420: }
1421: /* Copy constraint indices */
1422: for (subp = 0; subp < numSubpoints; ++subp) {
1423: PetscInt cdof;
1425: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1426: if (cdof) {
1427: for (f = 0; f < numFields; ++f) {
1428: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1429: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1430: }
1431: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1432: PetscSectionSetConstraintIndices(*subs, subp, indices);
1433: }
1434: }
1435: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1436: return(0);
1437: }
1439: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1440: {
1441: PetscInt p;
1442: PetscMPIInt rank;
1446: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1447: PetscViewerASCIIPushSynchronized(viewer);
1448: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1449: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1450: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1451: PetscInt b;
1453: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1454: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1455: PetscViewerASCIISynchronizedPrintf(viewer, " %d", s->bcIndices[s->bc->atlasOff[p]+b]);
1456: }
1457: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1458: } else {
1459: PetscViewerASCIISynchronizedPrintf(viewer, " (%4d) dim %2d offset %3d\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1460: }
1461: }
1462: PetscViewerFlush(viewer);
1463: PetscViewerASCIIPopSynchronized(viewer);
1464: if (s->sym) {
1465: PetscViewerASCIIPushTab(viewer);
1466: PetscSectionSymView(s->sym,viewer);
1467: PetscViewerASCIIPopTab(viewer);
1468: }
1469: return(0);
1470: }
1472: /*@C
1473: PetscSectionView - Views a PetscSection
1475: Collective on PetscSection
1477: Input Parameters:
1478: + s - the PetscSection object to view
1479: - v - the viewer
1481: Level: developer
1483: .seealso PetscSectionCreate(), PetscSectionDestroy()
1484: @*/
1485: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1486: {
1487: PetscBool isascii;
1488: PetscInt f;
1493: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1495: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1496: if (isascii) {
1497: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1498: if (s->numFields) {
1499: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1500: for (f = 0; f < s->numFields; ++f) {
1501: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
1502: PetscSectionView_ASCII(s->field[f], viewer);
1503: }
1504: } else {
1505: PetscSectionView_ASCII(s, viewer);
1506: }
1507: }
1508: return(0);
1509: }
1511: /*@
1512: PetscSectionReset - Frees all section data.
1514: Not collective
1516: Input Parameters:
1517: . s - the PetscSection
1519: Level: developer
1521: .seealso: PetscSection, PetscSectionCreate()
1522: @*/
1523: PetscErrorCode PetscSectionReset(PetscSection s)
1524: {
1525: PetscInt f;
1529: PetscFree(s->numFieldComponents);
1530: for (f = 0; f < s->numFields; ++f) {
1531: PetscSectionDestroy(&s->field[f]);
1532: PetscFree(s->fieldNames[f]);
1533: }
1534: PetscFree(s->fieldNames);
1535: PetscFree(s->field);
1536: PetscSectionDestroy(&s->bc);
1537: PetscFree(s->bcIndices);
1538: PetscFree2(s->atlasDof, s->atlasOff);
1539: PetscSectionDestroy(&s->clSection);
1540: ISDestroy(&s->clPoints);
1541: ISDestroy(&s->perm);
1542: PetscFree(s->clPerm);
1543: PetscFree(s->clInvPerm);
1544: PetscSectionSymDestroy(&s->sym);
1546: s->pStart = -1;
1547: s->pEnd = -1;
1548: s->maxDof = 0;
1549: s->setup = PETSC_FALSE;
1550: s->numFields = 0;
1551: s->clObj = NULL;
1552: return(0);
1553: }
1555: /*@
1556: PetscSectionDestroy - Frees a section object and frees its range if that exists.
1558: Not collective
1560: Input Parameters:
1561: . s - the PetscSection
1563: Level: developer
1565: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1566: recommended they not be used in user codes unless you really gain something in their use.
1568: .seealso: PetscSection, PetscSectionCreate()
1569: @*/
1570: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1571: {
1575: if (!*s) return(0);
1577: if (--((PetscObject)(*s))->refct > 0) {
1578: *s = NULL;
1579: return(0);
1580: }
1581: PetscSectionReset(*s);
1582: PetscHeaderDestroy(s);
1583: return(0);
1584: }
1586: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1587: {
1588: const PetscInt p = point - s->pStart;
1591: *values = &baseArray[s->atlasOff[p]];
1592: return(0);
1593: }
1595: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1596: {
1597: PetscInt *array;
1598: const PetscInt p = point - s->pStart;
1599: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1600: PetscInt cDim = 0;
1604: PetscSectionGetConstraintDof(s, p, &cDim);
1605: array = &baseArray[s->atlasOff[p]];
1606: if (!cDim) {
1607: if (orientation >= 0) {
1608: const PetscInt dim = s->atlasDof[p];
1609: PetscInt i;
1611: if (mode == INSERT_VALUES) {
1612: for (i = 0; i < dim; ++i) array[i] = values[i];
1613: } else {
1614: for (i = 0; i < dim; ++i) array[i] += values[i];
1615: }
1616: } else {
1617: PetscInt offset = 0;
1618: PetscInt j = -1, field, i;
1620: for (field = 0; field < s->numFields; ++field) {
1621: const PetscInt dim = s->field[field]->atlasDof[p];
1623: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1624: offset += dim;
1625: }
1626: }
1627: } else {
1628: if (orientation >= 0) {
1629: const PetscInt dim = s->atlasDof[p];
1630: PetscInt cInd = 0, i;
1631: const PetscInt *cDof;
1633: PetscSectionGetConstraintIndices(s, point, &cDof);
1634: if (mode == INSERT_VALUES) {
1635: for (i = 0; i < dim; ++i) {
1636: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1637: array[i] = values[i];
1638: }
1639: } else {
1640: for (i = 0; i < dim; ++i) {
1641: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
1642: array[i] += values[i];
1643: }
1644: }
1645: } else {
1646: const PetscInt *cDof;
1647: PetscInt offset = 0;
1648: PetscInt cOffset = 0;
1649: PetscInt j = 0, field;
1651: PetscSectionGetConstraintIndices(s, point, &cDof);
1652: for (field = 0; field < s->numFields; ++field) {
1653: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
1654: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
1655: const PetscInt sDim = dim - tDim;
1656: PetscInt cInd = 0, i,k;
1658: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
1659: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
1660: array[j] = values[k];
1661: }
1662: offset += dim;
1663: cOffset += dim - tDim;
1664: }
1665: }
1666: }
1667: return(0);
1668: }
1670: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
1671: {
1675: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
1676: return(0);
1677: }
1679: /*@C
1680: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
1682: Input Parameters:
1683: + s - The PetscSection
1684: - point - The point
1686: Output Parameter:
1687: . indices - The constrained dofs
1689: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
1691: Level: advanced
1693: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1694: @*/
1695: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
1696: {
1700: if (s->bc) {
1701: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
1702: } else *indices = NULL;
1703: return(0);
1704: }
1706: /*@C
1707: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
1709: Input Parameters:
1710: + s - The PetscSection
1711: . point - The point
1712: - indices - The constrained dofs
1714: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
1716: Level: advanced
1718: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
1719: @*/
1720: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
1721: {
1725: if (s->bc) {
1726: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
1727: }
1728: return(0);
1729: }
1731: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
1732: {
1736: 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);
1737: PetscSectionGetConstraintIndices(s->field[field], point, indices);
1738: return(0);
1739: }
1741: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
1742: {
1746: 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);
1747: PetscSectionSetConstraintIndices(s->field[field], point, indices);
1748: return(0);
1749: }
1751: /*@
1752: PetscSectionPermute - Reorder the section according to the input point permutation
1754: Collective on PetscSection
1756: Input Parameter:
1757: + section - The PetscSection object
1758: - perm - The point permutation, old point p becomes new point perm[p]
1760: Output Parameter:
1761: . sectionNew - The permuted PetscSection
1763: Level: intermediate
1765: .keywords: mesh
1766: .seealso: MatPermute()
1767: @*/
1768: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
1769: {
1770: PetscSection s = section, sNew;
1771: const PetscInt *perm;
1772: PetscInt numFields, f, numPoints, pStart, pEnd, p;
1773: PetscErrorCode ierr;
1779: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
1780: PetscSectionGetNumFields(s, &numFields);
1781: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
1782: for (f = 0; f < numFields; ++f) {
1783: const char *name;
1784: PetscInt numComp;
1786: PetscSectionGetFieldName(s, f, &name);
1787: PetscSectionSetFieldName(sNew, f, name);
1788: PetscSectionGetFieldComponents(s, f, &numComp);
1789: PetscSectionSetFieldComponents(sNew, f, numComp);
1790: }
1791: ISGetLocalSize(permutation, &numPoints);
1792: ISGetIndices(permutation, &perm);
1793: PetscSectionGetChart(s, &pStart, &pEnd);
1794: PetscSectionSetChart(sNew, pStart, pEnd);
1795: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
1796: for (p = pStart; p < pEnd; ++p) {
1797: PetscInt dof, cdof;
1799: PetscSectionGetDof(s, p, &dof);
1800: PetscSectionSetDof(sNew, perm[p], dof);
1801: PetscSectionGetConstraintDof(s, p, &cdof);
1802: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
1803: for (f = 0; f < numFields; ++f) {
1804: PetscSectionGetFieldDof(s, p, f, &dof);
1805: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
1806: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1807: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
1808: }
1809: }
1810: PetscSectionSetUp(sNew);
1811: for (p = pStart; p < pEnd; ++p) {
1812: const PetscInt *cind;
1813: PetscInt cdof;
1815: PetscSectionGetConstraintDof(s, p, &cdof);
1816: if (cdof) {
1817: PetscSectionGetConstraintIndices(s, p, &cind);
1818: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
1819: }
1820: for (f = 0; f < numFields; ++f) {
1821: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1822: if (cdof) {
1823: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
1824: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
1825: }
1826: }
1827: }
1828: ISRestoreIndices(permutation, &perm);
1829: *sectionNew = sNew;
1830: return(0);
1831: }
1833: /*@C
1834: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
1836: Collective
1838: Input Parameters:
1839: + sf - The SF
1840: - rootSection - Section defined on root space
1842: Output Parameters:
1843: + remoteOffsets - root offsets in leaf storage, or NULL
1844: - leafSection - Section defined on the leaf space
1846: Level: intermediate
1848: .seealso: PetscSFCreate()
1849: @*/
1850: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
1851: {
1852: PetscSF embedSF;
1853: const PetscInt *ilocal, *indices;
1854: IS selected;
1855: PetscInt numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;
1859: PetscSectionGetNumFields(rootSection, &numFields);
1860: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
1861: for (f = 0; f < numFields; ++f) {
1862: PetscInt numComp = 0;
1863: PetscSectionGetFieldComponents(rootSection, f, &numComp);
1864: PetscSectionSetFieldComponents(leafSection, f, numComp);
1865: }
1866: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1867: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
1868: rpEnd = PetscMin(rpEnd,nroots);
1869: rpEnd = PetscMax(rpStart,rpEnd);
1870: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1871: ISGetIndices(selected, &indices);
1872: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1873: ISRestoreIndices(selected, &indices);
1874: ISDestroy(&selected);
1875: PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
1876: if (nleaves && ilocal) {
1877: for (i = 0; i < nleaves; ++i) {
1878: lpStart = PetscMin(lpStart, ilocal[i]);
1879: lpEnd = PetscMax(lpEnd, ilocal[i]);
1880: }
1881: ++lpEnd;
1882: } else {
1883: lpStart = 0;
1884: lpEnd = nleaves;
1885: }
1886: PetscSectionSetChart(leafSection, lpStart, lpEnd);
1887: /* Could fuse these at the cost of a copy and extra allocation */
1888: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1889: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
1890: for (f = 0; f < numFields; ++f) {
1891: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1892: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
1893: }
1894: if (remoteOffsets) {
1895: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
1896: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1897: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1898: }
1899: PetscSFDestroy(&embedSF);
1900: PetscSectionSetUp(leafSection);
1901: return(0);
1902: }
1904: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
1905: {
1906: PetscSF embedSF;
1907: const PetscInt *indices;
1908: IS selected;
1909: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
1910: PetscErrorCode ierr;
1913: *remoteOffsets = NULL;
1914: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
1915: if (numRoots < 0) return(0);
1916: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
1917: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1918: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
1919: ISGetIndices(selected, &indices);
1920: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
1921: ISRestoreIndices(selected, &indices);
1922: ISDestroy(&selected);
1923: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
1924: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1925: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
1926: PetscSFDestroy(&embedSF);
1927: return(0);
1928: }
1930: /*@C
1931: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
1933: Input Parameters:
1934: + sf - The SF
1935: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
1936: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
1937: - leafSection - Data layout of local points for incoming data (this is the distributed section)
1939: Output Parameters:
1940: - sectionSF - The new SF
1942: Note: Either rootSection or remoteOffsets can be specified
1944: Level: intermediate
1946: .seealso: PetscSFCreate()
1947: @*/
1948: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1949: {
1950: MPI_Comm comm;
1951: const PetscInt *localPoints;
1952: const PetscSFNode *remotePoints;
1953: PetscInt lpStart, lpEnd;
1954: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
1955: PetscInt *localIndices;
1956: PetscSFNode *remoteIndices;
1957: PetscInt i, ind;
1958: PetscErrorCode ierr;
1966: PetscObjectGetComm((PetscObject)sf,&comm);
1967: PetscSFCreate(comm, sectionSF);
1968: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
1969: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
1970: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
1971: if (numRoots < 0) return(0);
1972: for (i = 0; i < numPoints; ++i) {
1973: PetscInt localPoint = localPoints ? localPoints[i] : i;
1974: PetscInt dof;
1976: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
1977: PetscSectionGetDof(leafSection, localPoint, &dof);
1978: numIndices += dof;
1979: }
1980: }
1981: PetscMalloc1(numIndices, &localIndices);
1982: PetscMalloc1(numIndices, &remoteIndices);
1983: /* Create new index graph */
1984: for (i = 0, ind = 0; i < numPoints; ++i) {
1985: PetscInt localPoint = localPoints ? localPoints[i] : i;
1986: PetscInt rank = remotePoints[i].rank;
1988: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
1989: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
1990: PetscInt loff, dof, d;
1992: PetscSectionGetOffset(leafSection, localPoint, &loff);
1993: PetscSectionGetDof(leafSection, localPoint, &dof);
1994: for (d = 0; d < dof; ++d, ++ind) {
1995: localIndices[ind] = loff+d;
1996: remoteIndices[ind].rank = rank;
1997: remoteIndices[ind].index = remoteOffset+d;
1998: }
1999: }
2000: }
2001: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %d should be %d", ind, numIndices);
2002: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2003: return(0);
2004: }
2006: /*@
2007: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2009: Input Parameters:
2010: + section - The PetscSection
2011: . obj - A PetscObject which serves as the key for this index
2012: . clSection - Section giving the size of the closure of each point
2013: - clPoints - IS giving the points in each closure
2015: Note: We compress out closure points with no dofs in this section
2017: Level: intermediate
2019: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2020: @*/
2021: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2022: {
2026: if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2027: section->clObj = obj;
2028: PetscSectionDestroy(§ion->clSection);
2029: ISDestroy(§ion->clPoints);
2030: section->clSection = clSection;
2031: section->clPoints = clPoints;
2032: return(0);
2033: }
2035: /*@
2036: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2038: Input Parameters:
2039: + section - The PetscSection
2040: - obj - A PetscObject which serves as the key for this index
2042: Output Parameters:
2043: + clSection - Section giving the size of the closure of each point
2044: - clPoints - IS giving the points in each closure
2046: Note: We compress out closure points with no dofs in this section
2048: Level: intermediate
2050: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2051: @*/
2052: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2053: {
2055: if (section->clObj == obj) {
2056: if (clSection) *clSection = section->clSection;
2057: if (clPoints) *clPoints = section->clPoints;
2058: } else {
2059: if (clSection) *clSection = NULL;
2060: if (clPoints) *clPoints = NULL;
2061: }
2062: return(0);
2063: }
2065: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2066: {
2067: PetscInt i;
2071: if (section->clObj != obj) {
2072: PetscSectionDestroy(§ion->clSection);
2073: ISDestroy(§ion->clPoints);
2074: }
2075: section->clObj = obj;
2076: PetscFree(section->clPerm);
2077: PetscFree(section->clInvPerm);
2078: section->clSize = clSize;
2079: if (mode == PETSC_COPY_VALUES) {
2080: PetscMalloc1(clSize, §ion->clPerm);
2081: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2082: PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2083: } else if (mode == PETSC_OWN_POINTER) {
2084: section->clPerm = clPerm;
2085: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2086: PetscMalloc1(clSize, §ion->clInvPerm);
2087: for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2088: return(0);
2089: }
2091: /*@
2092: PetscSectionSetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2094: Input Parameters:
2095: + section - The PetscSection
2096: . obj - A PetscObject which serves as the key for this index
2097: - perm - Permutation of the cell dof closure
2099: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2100: other points (like faces).
2102: Level: intermediate
2104: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2105: @*/
2106: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2107: {
2108: const PetscInt *clPerm = NULL;
2109: PetscInt clSize = 0;
2110: PetscErrorCode ierr;
2113: if (perm) {
2114: ISGetLocalSize(perm, &clSize);
2115: ISGetIndices(perm, &clPerm);
2116: }
2117: PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2118: if (perm) {ISRestoreIndices(perm, &clPerm);}
2119: return(0);
2120: }
2122: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2123: {
2125: if (section->clObj == obj) {
2126: if (size) *size = section->clSize;
2127: if (perm) *perm = section->clPerm;
2128: } else {
2129: if (size) *size = 0;
2130: if (perm) *perm = NULL;
2131: }
2132: return(0);
2133: }
2135: /*@
2136: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2138: Input Parameters:
2139: + section - The PetscSection
2140: - obj - A PetscObject which serves as the key for this index
2142: Output Parameter:
2143: . perm - The dof closure permutation
2145: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2146: other points (like faces).
2148: The user must destroy the IS that is returned.
2150: Level: intermediate
2152: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2153: @*/
2154: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2155: {
2156: const PetscInt *clPerm;
2157: PetscInt clSize;
2158: PetscErrorCode ierr;
2161: PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2162: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2163: return(0);
2164: }
2166: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2167: {
2169: if (section->clObj == obj) {
2170: if (size) *size = section->clSize;
2171: if (perm) *perm = section->clInvPerm;
2172: } else {
2173: if (size) *size = 0;
2174: if (perm) *perm = NULL;
2175: }
2176: return(0);
2177: }
2179: /*@
2180: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2182: Input Parameters:
2183: + section - The PetscSection
2184: - obj - A PetscObject which serves as the key for this index
2186: Output Parameters:
2187: + size - The dof closure size
2188: - perm - The dof closure permutation
2190: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2191: other points (like faces).
2193: The user must destroy the IS that is returned.
2195: Level: intermediate
2197: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2198: @*/
2199: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2200: {
2201: const PetscInt *clPerm;
2202: PetscInt clSize;
2203: PetscErrorCode ierr;
2206: PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2207: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2208: return(0);
2209: }
2211: /*@
2212: PetscSectionGetField - Get the subsection associated with a single field
2214: Input Parameters:
2215: + s - The PetscSection
2216: - field - The field number
2218: Output Parameter:
2219: . subs - The subsection for the given field
2221: Level: intermediate
2223: .seealso: PetscSectionSetNumFields()
2224: @*/
2225: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2226: {
2232: 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);
2233: PetscObjectReference((PetscObject) s->field[field]);
2234: *subs = s->field[field];
2235: return(0);
2236: }
2238: PetscClassId PETSC_SECTION_SYM_CLASSID;
2239: PetscFunctionList PetscSectionSymList = NULL;
2241: /*@
2242: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
2244: Collective on MPI_Comm
2246: Input Parameter:
2247: . comm - the MPI communicator
2249: Output Parameter:
2250: . sym - pointer to the new set of symmetries
2252: Level: developer
2254: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2255: @*/
2256: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2257: {
2262: ISInitializePackage();
2263: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2264: return(0);
2265: }
2267: /*@C
2268: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
2270: Collective on PetscSectionSym
2272: Input Parameters:
2273: + sym - The section symmetry object
2274: - method - The name of the section symmetry type
2276: Level: developer
2278: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2279: @*/
2280: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2281: {
2282: PetscErrorCode (*r)(PetscSectionSym);
2283: PetscBool match;
2288: PetscObjectTypeCompare((PetscObject) sym, method, &match);
2289: if (match) return(0);
2291: PetscFunctionListFind(PetscSectionSymList,method,&r);
2292: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2293: if (sym->ops->destroy) {
2294: (*sym->ops->destroy)(sym);
2295: sym->ops->destroy = NULL;
2296: }
2297: (*r)(sym);
2298: PetscObjectChangeTypeName((PetscObject)sym,method);
2299: return(0);
2300: }
2303: /*@C
2304: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
2306: Not Collective
2308: Input Parameter:
2309: . sym - The section symmetry
2311: Output Parameter:
2312: . type - The index set type name
2314: Level: developer
2316: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2317: @*/
2318: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2319: {
2323: *type = ((PetscObject)sym)->type_name;
2324: return(0);
2325: }
2327: /*@C
2328: PetscSectionSymRegister - Adds a new section symmetry implementation
2330: Not Collective
2332: Input Parameters:
2333: + name - The name of a new user-defined creation routine
2334: - create_func - The creation routine itself
2336: Notes:
2337: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
2339: Level: developer
2341: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2342: @*/
2343: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2344: {
2348: ISInitializePackage();
2349: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2350: return(0);
2351: }
2353: /*@
2354: PetscSectionSymDestroy - Destroys a section symmetry.
2356: Collective on PetscSectionSym
2358: Input Parameters:
2359: . sym - the section symmetry
2361: Level: developer
2363: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2364: @*/
2365: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2366: {
2367: SymWorkLink link,next;
2371: if (!*sym) return(0);
2373: if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2374: if ((*sym)->ops->destroy) {
2375: (*(*sym)->ops->destroy)(*sym);
2376: }
2377: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2378: for (link=(*sym)->workin; link; link=next) {
2379: next = link->next;
2380: PetscFree2(link->perms,link->rots);
2381: PetscFree(link);
2382: }
2383: (*sym)->workin = NULL;
2384: PetscHeaderDestroy(sym);
2385: return(0);
2386: }
2388: /*@C
2389: PetscSectionSymView - Displays a section symmetry
2391: Collective on PetscSectionSym
2393: Input Parameters:
2394: + sym - the index set
2395: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
2397: Level: developer
2399: .seealso: PetscViewerASCIIOpen()
2400: @*/
2401: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2402: {
2407: if (!viewer) {
2408: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2409: }
2412: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2413: if (sym->ops->view) {
2414: (*sym->ops->view)(sym,viewer);
2415: }
2416: return(0);
2417: }
2419: /*@
2420: PetscSectionSetSym - Set the symmetries for the data referred to by the section
2422: Collective on PetscSection
2424: Input Parameters:
2425: + section - the section describing data layout
2426: - sym - the symmetry describing the affect of orientation on the access of the data
2428: Level: developer
2430: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2431: @*/
2432: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2433: {
2440: PetscObjectReference((PetscObject)sym);
2441: PetscSectionSymDestroy(&(section->sym));
2442: section->sym = sym;
2443: return(0);
2444: }
2446: /*@
2447: PetscSectionGetSym - Get the symmetries for the data referred to by the section
2449: Not collective
2451: Input Parameters:
2452: . section - the section describing data layout
2454: Output Parameters:
2455: . sym - the symmetry describing the affect of orientation on the access of the data
2457: Level: developer
2459: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2460: @*/
2461: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2462: {
2465: *sym = section->sym;
2466: return(0);
2467: }
2469: /*@
2470: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
2472: Collective on PetscSection
2474: Input Parameters:
2475: + section - the section describing data layout
2476: . field - the field number
2477: - sym - the symmetry describing the affect of orientation on the access of the data
2479: Level: developer
2481: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2482: @*/
2483: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2484: {
2489: 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);
2490: PetscSectionSetSym(section->field[field],sym);
2491: return(0);
2492: }
2494: /*@
2495: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
2497: Collective on PetscSection
2499: Input Parameters:
2500: + section - the section describing data layout
2501: - field - the field number
2503: Output Parameters:
2504: . sym - the symmetry describing the affect of orientation on the access of the data
2506: Level: developer
2508: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2509: @*/
2510: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2511: {
2514: 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);
2515: *sym = section->field[field]->sym;
2516: return(0);
2517: }
2519: /*@C
2520: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
2522: Not collective
2524: Input Parameters:
2525: + section - the section
2526: . numPoints - the number of points
2527: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2528: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2529: context, see DMPlexGetConeOrientation()).
2531: Output Parameter:
2532: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2533: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2534: identity).
2536: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2537: .vb
2538: const PetscInt **perms;
2539: const PetscScalar **rots;
2540: PetscInt lOffset;
2542: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2543: for (i = 0, lOffset = 0; i < numPoints; i++) {
2544: PetscInt point = points[2*i], dof, sOffset;
2545: const PetscInt *perm = perms ? perms[i] : NULL;
2546: const PetscScalar *rot = rots ? rots[i] : NULL;
2548: PetscSectionGetDof(section,point,&dof);
2549: PetscSectionGetOffset(section,point,&sOffset);
2551: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
2552: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
2553: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
2554: lOffset += dof;
2555: }
2556: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2557: .ve
2559: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2560: .vb
2561: const PetscInt **perms;
2562: const PetscScalar **rots;
2563: PetscInt lOffset;
2565: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2566: for (i = 0, lOffset = 0; i < numPoints; i++) {
2567: PetscInt point = points[2*i], dof, sOffset;
2568: const PetscInt *perm = perms ? perms[i] : NULL;
2569: const PetscScalar *rot = rots ? rots[i] : NULL;
2571: PetscSectionGetDof(section,point,&dof);
2572: PetscSectionGetOffset(section,point,&sOff);
2574: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2575: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
2576: offset += dof;
2577: }
2578: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2579: .ve
2581: Level: developer
2583: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2584: @*/
2585: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2586: {
2587: PetscSectionSym sym;
2588: PetscErrorCode ierr;
2593: if (perms) *perms = NULL;
2594: if (rots) *rots = NULL;
2595: sym = section->sym;
2596: if (sym && (perms || rots)) {
2597: SymWorkLink link;
2599: if (sym->workin) {
2600: link = sym->workin;
2601: sym->workin = sym->workin->next;
2602: } else {
2603: PetscNewLog(sym,&link);
2604: }
2605: if (numPoints > link->numPoints) {
2606: PetscFree2(link->perms,link->rots);
2607: PetscMalloc2(numPoints,&link->perms,numPoints,&link->rots);
2608: link->numPoints = numPoints;
2609: }
2610: link->next = sym->workout;
2611: sym->workout = link;
2612: PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2613: PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2614: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2615: if (perms) *perms = link->perms;
2616: if (rots) *rots = link->rots;
2617: }
2618: return(0);
2619: }
2621: /*@C
2622: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
2624: Not collective
2626: Input Parameters:
2627: + section - the section
2628: . numPoints - the number of points
2629: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2630: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2631: context, see DMPlexGetConeOrientation()).
2633: Output Parameter:
2634: + perms - The permutations for the given orientations: set to NULL at conclusion
2635: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
2637: Level: developer
2639: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2640: @*/
2641: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2642: {
2643: PetscSectionSym sym;
2647: sym = section->sym;
2648: if (sym && (perms || rots)) {
2649: SymWorkLink *p,link;
2651: for (p=&sym->workout; (link=*p); p=&link->next) {
2652: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
2653: *p = link->next;
2654: link->next = sym->workin;
2655: sym->workin = link;
2656: if (perms) *perms = NULL;
2657: if (rots) *rots = NULL;
2658: return(0);
2659: }
2660: }
2661: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
2662: }
2663: return(0);
2664: }
2666: /*@C
2667: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
2669: Not collective
2671: Input Parameters:
2672: + section - the section
2673: . field - the field of the section
2674: . numPoints - the number of points
2675: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2676: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2677: context, see DMPlexGetConeOrientation()).
2679: Output Parameter:
2680: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2681: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2682: identity).
2684: Level: developer
2686: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
2687: @*/
2688: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2689: {
2694: 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);
2695: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
2696: return(0);
2697: }
2699: /*@C
2700: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
2702: Not collective
2704: Input Parameters:
2705: + section - the section
2706: . field - the field number
2707: . numPoints - the number of points
2708: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2709: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2710: context, see DMPlexGetConeOrientation()).
2712: Output Parameter:
2713: + perms - The permutations for the given orientations: set to NULL at conclusion
2714: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
2716: Level: developer
2718: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2719: @*/
2720: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2721: {
2726: 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);
2727: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
2728: return(0);
2729: }