Actual source code: section.c
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsection.h>
7: #include <petscsf.h>
8: #include <petscviewer.h>
10: PetscClassId PETSC_SECTION_CLASSID;
12: /*@
13: PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
15: Collective
17: Input Parameters:
18: + comm - the MPI communicator
19: - s - pointer to the section
21: Level: beginner
23: Notes:
24: Typical calling sequence
25: $ PetscSectionCreate(MPI_Comm,PetscSection *);
26: $ PetscSectionSetNumFields(PetscSection, numFields);
27: $ PetscSectionSetChart(PetscSection,low,high);
28: $ PetscSectionSetDof(PetscSection,point,numdof);
29: $ PetscSectionSetUp(PetscSection);
30: $ PetscSectionGetOffset(PetscSection,point,PetscInt *);
31: $ PetscSectionDestroy(PetscSection);
33: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
34: recommended they not be used in user codes unless you really gain something in their use.
36: .seealso: PetscSection, PetscSectionDestroy()
37: @*/
38: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
39: {
44: ISInitializePackage();
46: PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
48: (*s)->pStart = -1;
49: (*s)->pEnd = -1;
50: (*s)->perm = NULL;
51: (*s)->pointMajor = PETSC_TRUE;
52: (*s)->maxDof = 0;
53: (*s)->atlasDof = NULL;
54: (*s)->atlasOff = NULL;
55: (*s)->bc = NULL;
56: (*s)->bcIndices = NULL;
57: (*s)->setup = PETSC_FALSE;
58: (*s)->numFields = 0;
59: (*s)->fieldNames = NULL;
60: (*s)->field = NULL;
61: (*s)->useFieldOff = PETSC_FALSE;
62: (*s)->compNames = NULL;
63: (*s)->clObj = NULL;
64: (*s)->clHash = NULL;
65: (*s)->clSection = NULL;
66: (*s)->clPoints = NULL;
67: return(0);
68: }
70: /*@
71: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
73: Collective
75: Input Parameter:
76: . section - the PetscSection
78: Output Parameter:
79: . newSection - the copy
81: Level: intermediate
83: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
84: @*/
85: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
86: {
87: PetscSectionSym sym;
88: IS perm;
89: PetscInt numFields, f, c, pStart, pEnd, p;
90: PetscErrorCode ierr;
95: PetscSectionReset(newSection);
96: PetscSectionGetNumFields(section, &numFields);
97: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
98: for (f = 0; f < numFields; ++f) {
99: const char *fieldName = NULL, *compName = NULL;
100: PetscInt numComp = 0;
102: PetscSectionGetFieldName(section, f, &fieldName);
103: PetscSectionSetFieldName(newSection, f, fieldName);
104: PetscSectionGetFieldComponents(section, f, &numComp);
105: PetscSectionSetFieldComponents(newSection, f, numComp);
106: for (c = 0; c < numComp; ++c) {
107: PetscSectionGetComponentName(section, f, c, &compName);
108: PetscSectionSetComponentName(newSection, f, c, compName);
109: }
110: PetscSectionGetFieldSym(section, f, &sym);
111: PetscSectionSetFieldSym(newSection, f, sym);
112: }
113: PetscSectionGetChart(section, &pStart, &pEnd);
114: PetscSectionSetChart(newSection, pStart, pEnd);
115: PetscSectionGetPermutation(section, &perm);
116: PetscSectionSetPermutation(newSection, perm);
117: PetscSectionGetSym(section, &sym);
118: PetscSectionSetSym(newSection, sym);
119: for (p = pStart; p < pEnd; ++p) {
120: PetscInt dof, cdof, fcdof = 0;
122: PetscSectionGetDof(section, p, &dof);
123: PetscSectionSetDof(newSection, p, dof);
124: PetscSectionGetConstraintDof(section, p, &cdof);
125: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
126: for (f = 0; f < numFields; ++f) {
127: PetscSectionGetFieldDof(section, p, f, &dof);
128: PetscSectionSetFieldDof(newSection, p, f, dof);
129: if (cdof) {
130: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
131: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
132: }
133: }
134: }
135: PetscSectionSetUp(newSection);
136: for (p = pStart; p < pEnd; ++p) {
137: PetscInt off, cdof, fcdof = 0;
138: const PetscInt *cInd;
140: /* Must set offsets in case they do not agree with the prefix sums */
141: PetscSectionGetOffset(section, p, &off);
142: PetscSectionSetOffset(newSection, p, off);
143: PetscSectionGetConstraintDof(section, p, &cdof);
144: if (cdof) {
145: PetscSectionGetConstraintIndices(section, p, &cInd);
146: PetscSectionSetConstraintIndices(newSection, p, cInd);
147: for (f = 0; f < numFields; ++f) {
148: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
149: if (fcdof) {
150: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
151: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
152: }
153: }
154: }
155: }
156: return(0);
157: }
159: /*@
160: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
162: Collective
164: Input Parameter:
165: . section - the PetscSection
167: Output Parameter:
168: . newSection - the copy
170: Level: beginner
172: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
173: @*/
174: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
175: {
181: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
182: PetscSectionCopy(section, *newSection);
183: return(0);
184: }
186: /*@
187: PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database
189: Collective on PetscSection
191: Input Parameter:
192: . section - the PetscSection
194: Options Database:
195: . -petscsection_point_major the dof order
197: Level: intermediate
199: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
200: @*/
201: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
202: {
207: PetscObjectOptionsBegin((PetscObject) s);
208: PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
209: /* process any options handlers added with PetscObjectAddOptionsHandler() */
210: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
211: PetscOptionsEnd();
212: PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
213: return(0);
214: }
216: /*@
217: PetscSectionCompare - Compares two sections
219: Collective on PetscSection
221: Input Parameters:
222: + s1 - the first PetscSection
223: - s2 - the second PetscSection
225: Output Parameter:
226: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
228: Level: intermediate
230: Notes:
231: Field names are disregarded.
233: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
234: @*/
235: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
236: {
237: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
238: const PetscInt *idx1, *idx2;
239: IS perm1, perm2;
240: PetscBool flg;
241: PetscMPIInt mflg;
248: flg = PETSC_FALSE;
250: MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
251: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
252: *congruent = PETSC_FALSE;
253: return(0);
254: }
256: PetscSectionGetChart(s1, &pStart, &pEnd);
257: PetscSectionGetChart(s2, &n1, &n2);
258: if (pStart != n1 || pEnd != n2) goto not_congruent;
260: PetscSectionGetPermutation(s1, &perm1);
261: PetscSectionGetPermutation(s2, &perm2);
262: if (perm1 && perm2) {
263: ISEqual(perm1, perm2, congruent);
264: if (!(*congruent)) goto not_congruent;
265: } else if (perm1 != perm2) goto not_congruent;
267: for (p = pStart; p < pEnd; ++p) {
268: PetscSectionGetOffset(s1, p, &n1);
269: PetscSectionGetOffset(s2, p, &n2);
270: if (n1 != n2) goto not_congruent;
272: PetscSectionGetDof(s1, p, &n1);
273: PetscSectionGetDof(s2, p, &n2);
274: if (n1 != n2) goto not_congruent;
276: PetscSectionGetConstraintDof(s1, p, &ncdof);
277: PetscSectionGetConstraintDof(s2, p, &n2);
278: if (ncdof != n2) goto not_congruent;
280: PetscSectionGetConstraintIndices(s1, p, &idx1);
281: PetscSectionGetConstraintIndices(s2, p, &idx2);
282: PetscArraycmp(idx1, idx2, ncdof, congruent);
283: if (!(*congruent)) goto not_congruent;
284: }
286: PetscSectionGetNumFields(s1, &nfields);
287: PetscSectionGetNumFields(s2, &n2);
288: if (nfields != n2) goto not_congruent;
290: for (f = 0; f < nfields; ++f) {
291: PetscSectionGetFieldComponents(s1, f, &n1);
292: PetscSectionGetFieldComponents(s2, f, &n2);
293: if (n1 != n2) goto not_congruent;
295: for (p = pStart; p < pEnd; ++p) {
296: PetscSectionGetFieldOffset(s1, p, f, &n1);
297: PetscSectionGetFieldOffset(s2, p, f, &n2);
298: if (n1 != n2) goto not_congruent;
300: PetscSectionGetFieldDof(s1, p, f, &n1);
301: PetscSectionGetFieldDof(s2, p, f, &n2);
302: if (n1 != n2) goto not_congruent;
304: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
305: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
306: if (nfcdof != n2) goto not_congruent;
308: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
309: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
310: PetscArraycmp(idx1, idx2, nfcdof, congruent);
311: if (!(*congruent)) goto not_congruent;
312: }
313: }
315: flg = PETSC_TRUE;
316: not_congruent:
317: MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
318: return(0);
319: }
321: /*@
322: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
324: Not collective
326: Input Parameter:
327: . s - the PetscSection
329: Output Parameter:
330: . numFields - the number of fields defined, or 0 if none were defined
332: Level: intermediate
334: .seealso: PetscSectionSetNumFields()
335: @*/
336: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
337: {
341: *numFields = s->numFields;
342: return(0);
343: }
345: /*@
346: PetscSectionSetNumFields - Sets the number of fields.
348: Not collective
350: Input Parameters:
351: + s - the PetscSection
352: - numFields - the number of fields
354: Level: intermediate
356: .seealso: PetscSectionGetNumFields()
357: @*/
358: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
359: {
360: PetscInt f;
365: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %D must be positive", numFields);
366: PetscSectionReset(s);
368: s->numFields = numFields;
369: PetscMalloc1(s->numFields, &s->numFieldComponents);
370: PetscMalloc1(s->numFields, &s->fieldNames);
371: PetscMalloc1(s->numFields, &s->compNames);
372: PetscMalloc1(s->numFields, &s->field);
373: for (f = 0; f < s->numFields; ++f) {
374: char name[64];
376: s->numFieldComponents[f] = 1;
378: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
379: PetscSNPrintf(name, 64, "Field_%D", f);
380: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
381: PetscSNPrintf(name, 64, "Component_0");
382: PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
383: PetscStrallocpy(name, (char **) &s->compNames[f][0]);
384: }
385: return(0);
386: }
388: /*@C
389: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
391: Not Collective
393: Input Parameters:
394: + s - the PetscSection
395: - field - the field number
397: Output Parameter:
398: . fieldName - the field name
400: Level: intermediate
402: .seealso: PetscSectionSetFieldName()
403: @*/
404: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
405: {
409: 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);
410: *fieldName = s->fieldNames[field];
411: return(0);
412: }
414: /*@C
415: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
417: Not Collective
419: Input Parameters:
420: + s - the PetscSection
421: . field - the field number
422: - fieldName - the field name
424: Level: intermediate
426: .seealso: PetscSectionGetFieldName()
427: @*/
428: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
429: {
435: 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);
436: PetscFree(s->fieldNames[field]);
437: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
438: return(0);
439: }
441: /*@C
442: PetscSectionGetComponentName - Gets the name of a field component in the PetscSection
444: Not Collective
446: Input Parameters:
447: + s - the PetscSection
448: . field - the field number
449: . comp - the component number
450: - compName - the component name
452: Level: intermediate
454: .seealso: PetscSectionSetComponentName()
455: @*/
456: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
457: {
461: 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);
462: if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
463: *compName = s->compNames[field][comp];
464: return(0);
465: }
467: /*@C
468: PetscSectionSetComponentName - Sets the name of a field component in the PetscSection
470: Not Collective
472: Input Parameters:
473: + s - the PetscSection
474: . field - the field number
475: . comp - the component number
476: - compName - the component name
478: Level: intermediate
480: .seealso: PetscSectionGetComponentName()
481: @*/
482: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
483: {
489: 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);
490: if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
491: PetscFree(s->compNames[field][comp]);
492: PetscStrallocpy(compName, (char**) &s->compNames[field][comp]);
493: return(0);
494: }
496: /*@
497: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
499: Not collective
501: Input Parameters:
502: + s - the PetscSection
503: - field - the field number
505: Output Parameter:
506: . numComp - the number of field components
508: Level: intermediate
510: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
511: @*/
512: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
513: {
517: 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);
518: *numComp = s->numFieldComponents[field];
519: return(0);
520: }
522: /*@
523: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
525: Not collective
527: Input Parameters:
528: + s - the PetscSection
529: . field - the field number
530: - numComp - the number of field components
532: Level: intermediate
534: .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
535: @*/
536: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
537: {
539: PetscInt c;
540: char name[64];
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: if (s->compNames) {
546: for (c = 0; c < s->numFieldComponents[field]; ++c) {
547: PetscFree(s->compNames[field][c]);
548: }
549: PetscFree(s->compNames[field]);
550: }
552: s->numFieldComponents[field] = numComp;
553: if (numComp) {
554: PetscMalloc1(numComp, (char ***) &s->compNames[field]);
555: for (c = 0; c < numComp; ++c) {
556: PetscSNPrintf(name, 64, "%D", c);
557: PetscStrallocpy(name, (char **) &s->compNames[field][c]);
558: }
559: }
560: return(0);
561: }
563: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
564: {
568: if (!s->bc) {
569: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
570: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
571: }
572: return(0);
573: }
575: /*@
576: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
578: Not collective
580: Input Parameter:
581: . s - the PetscSection
583: Output Parameters:
584: + pStart - the first point
585: - pEnd - one past the last point
587: Level: intermediate
589: .seealso: PetscSectionSetChart(), PetscSectionCreate()
590: @*/
591: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
592: {
595: if (pStart) *pStart = s->pStart;
596: if (pEnd) *pEnd = s->pEnd;
597: return(0);
598: }
600: /*@
601: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
603: Not collective
605: Input Parameters:
606: + s - the PetscSection
607: . pStart - the first point
608: - pEnd - one past the last point
610: Level: intermediate
612: .seealso: PetscSectionGetChart(), PetscSectionCreate()
613: @*/
614: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
615: {
616: PetscInt f;
621: /* Cannot Reset() because it destroys field information */
622: s->setup = PETSC_FALSE;
623: PetscSectionDestroy(&s->bc);
624: PetscFree(s->bcIndices);
625: PetscFree2(s->atlasDof, s->atlasOff);
627: s->pStart = pStart;
628: s->pEnd = pEnd;
629: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
630: PetscArrayzero(s->atlasDof, pEnd - pStart);
631: for (f = 0; f < s->numFields; ++f) {
632: PetscSectionSetChart(s->field[f], pStart, pEnd);
633: }
634: return(0);
635: }
637: /*@
638: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
640: Not collective
642: Input Parameter:
643: . s - the PetscSection
645: Output Parameters:
646: . perm - The permutation as an IS
648: Level: intermediate
650: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
651: @*/
652: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
653: {
657: return(0);
658: }
660: /*@
661: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
663: Not collective
665: Input Parameters:
666: + s - the PetscSection
667: - perm - the permutation of points
669: Level: intermediate
671: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
672: @*/
673: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
674: {
680: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
681: if (s->perm != perm) {
682: ISDestroy(&s->perm);
683: if (perm) {
684: s->perm = perm;
685: PetscObjectReference((PetscObject) s->perm);
686: }
687: }
688: return(0);
689: }
691: /*@
692: PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major
694: Not collective
696: Input Parameter:
697: . s - the PetscSection
699: Output Parameter:
700: . pm - the flag for point major ordering
702: Level: intermediate
704: .seealso: PetscSectionSetPointMajor()
705: @*/
706: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
707: {
711: *pm = s->pointMajor;
712: return(0);
713: }
715: /*@
716: PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major
718: Not collective
720: Input Parameters:
721: + s - the PetscSection
722: - pm - the flag for point major ordering
724: Not collective
726: Level: intermediate
728: .seealso: PetscSectionGetPointMajor()
729: @*/
730: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
731: {
734: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
735: s->pointMajor = pm;
736: return(0);
737: }
739: /*@
740: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
742: Not collective
744: Input Parameters:
745: + s - the PetscSection
746: - point - the point
748: Output Parameter:
749: . numDof - the number of dof
751: Level: intermediate
753: .seealso: PetscSectionSetDof(), PetscSectionCreate()
754: @*/
755: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
756: {
760: if (PetscDefined(USE_DEBUG)) {
761: 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);
762: }
763: *numDof = s->atlasDof[point - s->pStart];
764: return(0);
765: }
767: /*@
768: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
770: Not collective
772: Input Parameters:
773: + s - the PetscSection
774: . point - the point
775: - numDof - the number of dof
777: Level: intermediate
779: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
780: @*/
781: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
782: {
785: 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);
786: s->atlasDof[point - s->pStart] = numDof;
787: return(0);
788: }
790: /*@
791: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
793: Not collective
795: Input Parameters:
796: + s - the PetscSection
797: . point - the point
798: - numDof - the number of additional dof
800: Level: intermediate
802: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
803: @*/
804: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
805: {
808: if (PetscDefined(USE_DEBUG)) {
809: 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);
810: }
811: s->atlasDof[point - s->pStart] += numDof;
812: return(0);
813: }
815: /*@
816: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given 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
828: Level: intermediate
830: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
831: @*/
832: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
833: {
838: 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);
839: PetscSectionGetDof(s->field[field], point, numDof);
840: return(0);
841: }
843: /*@
844: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
846: Not collective
848: Input Parameters:
849: + s - the PetscSection
850: . point - the point
851: . field - the field
852: - numDof - the number of dof
854: Level: intermediate
856: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
857: @*/
858: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
859: {
864: 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);
865: PetscSectionSetDof(s->field[field], point, numDof);
866: return(0);
867: }
869: /*@
870: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
872: Not collective
874: Input Parameters:
875: + s - the PetscSection
876: . point - the point
877: . field - the field
878: - numDof - the number of dof
880: Level: intermediate
882: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
883: @*/
884: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
885: {
890: 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);
891: PetscSectionAddDof(s->field[field], point, numDof);
892: return(0);
893: }
895: /*@
896: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
898: Not collective
900: Input Parameters:
901: + s - the PetscSection
902: - point - the point
904: Output Parameter:
905: . numDof - the number of dof which are fixed by constraints
907: Level: intermediate
909: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
910: @*/
911: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
912: {
918: if (s->bc) {
919: PetscSectionGetDof(s->bc, point, numDof);
920: } else *numDof = 0;
921: return(0);
922: }
924: /*@
925: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
927: Not collective
929: Input Parameters:
930: + s - the PetscSection
931: . point - the point
932: - numDof - the number of dof which are fixed by constraints
934: Level: intermediate
936: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
937: @*/
938: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
939: {
944: if (numDof) {
945: PetscSectionCheckConstraints_Static(s);
946: PetscSectionSetDof(s->bc, point, numDof);
947: }
948: return(0);
949: }
951: /*@
952: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
954: Not collective
956: Input Parameters:
957: + s - the PetscSection
958: . point - the point
959: - numDof - the number of additional dof which are fixed by constraints
961: Level: intermediate
963: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
964: @*/
965: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
966: {
971: if (numDof) {
972: PetscSectionCheckConstraints_Static(s);
973: PetscSectionAddDof(s->bc, point, numDof);
974: }
975: return(0);
976: }
978: /*@
979: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
981: Not collective
983: Input Parameters:
984: + s - the PetscSection
985: . point - the point
986: - field - the field
988: Output Parameter:
989: . numDof - the number of dof which are fixed by constraints
991: Level: intermediate
993: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
994: @*/
995: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
996: {
1002: 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);
1003: PetscSectionGetConstraintDof(s->field[field], point, numDof);
1004: return(0);
1005: }
1007: /*@
1008: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1010: Not collective
1012: Input Parameters:
1013: + s - the PetscSection
1014: . point - the point
1015: . field - the field
1016: - numDof - the number of dof which are fixed by constraints
1018: Level: intermediate
1020: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1021: @*/
1022: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1023: {
1028: 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);
1029: PetscSectionSetConstraintDof(s->field[field], point, numDof);
1030: return(0);
1031: }
1033: /*@
1034: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1036: Not collective
1038: Input Parameters:
1039: + s - the PetscSection
1040: . point - the point
1041: . field - the field
1042: - numDof - the number of additional dof which are fixed by constraints
1044: Level: intermediate
1046: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1047: @*/
1048: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1049: {
1054: 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);
1055: PetscSectionAddConstraintDof(s->field[field], point, numDof);
1056: return(0);
1057: }
1059: /*@
1060: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1062: Not collective
1064: Input Parameter:
1065: . s - the PetscSection
1067: Level: advanced
1069: .seealso: PetscSectionSetUp(), PetscSectionCreate()
1070: @*/
1071: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1072: {
1077: if (s->bc) {
1078: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
1080: PetscSectionSetUp(s->bc);
1081: PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1082: }
1083: return(0);
1084: }
1086: /*@
1087: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
1089: Not collective
1091: Input Parameter:
1092: . s - the PetscSection
1094: Level: intermediate
1096: .seealso: PetscSectionCreate()
1097: @*/
1098: PetscErrorCode PetscSectionSetUp(PetscSection s)
1099: {
1100: const PetscInt *pind = NULL;
1101: PetscInt offset = 0, foff, p, f;
1102: PetscErrorCode ierr;
1106: if (s->setup) return(0);
1107: s->setup = PETSC_TRUE;
1108: /* Set offsets and field offsets for all points */
1109: /* Assume that all fields have the same chart */
1110: if (s->perm) {ISGetIndices(s->perm, &pind);}
1111: if (s->pointMajor) {
1112: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1113: const PetscInt q = pind ? pind[p] : p;
1115: /* Set point offset */
1116: s->atlasOff[q] = offset;
1117: offset += s->atlasDof[q];
1118: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
1119: /* Set field offset */
1120: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1121: PetscSection sf = s->field[f];
1123: sf->atlasOff[q] = foff;
1124: foff += sf->atlasDof[q];
1125: }
1126: }
1127: } else {
1128: /* Set field offsets for all points */
1129: for (f = 0; f < s->numFields; ++f) {
1130: PetscSection sf = s->field[f];
1132: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1133: const PetscInt q = pind ? pind[p] : p;
1135: sf->atlasOff[q] = offset;
1136: offset += sf->atlasDof[q];
1137: }
1138: }
1139: /* Disable point offsets since these are unused */
1140: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1141: s->atlasOff[p] = -1;
1142: s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1143: }
1144: }
1145: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1146: /* Setup BC sections */
1147: PetscSectionSetUpBC(s);
1148: for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1149: return(0);
1150: }
1152: /*@
1153: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1155: Not collective
1157: Input Parameters:
1158: . s - the PetscSection
1160: Output Parameter:
1161: . maxDof - the maximum dof
1163: Level: intermediate
1165: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1166: @*/
1167: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1168: {
1172: *maxDof = s->maxDof;
1173: return(0);
1174: }
1176: /*@
1177: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
1179: Not collective
1181: Input Parameter:
1182: . s - the PetscSection
1184: Output Parameter:
1185: . size - the size of an array which can hold all the dofs
1187: Level: intermediate
1189: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1190: @*/
1191: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1192: {
1193: PetscInt p, n = 0;
1198: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1199: *size = n;
1200: return(0);
1201: }
1203: /*@
1204: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1206: Not collective
1208: Input Parameter:
1209: . s - the PetscSection
1211: Output Parameter:
1212: . size - the size of an array which can hold all unconstrained dofs
1214: Level: intermediate
1216: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1217: @*/
1218: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1219: {
1220: PetscInt p, n = 0;
1225: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1226: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1227: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1228: }
1229: *size = n;
1230: return(0);
1231: }
1233: /*@
1234: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1235: the local section and an SF describing the section point overlap.
1237: Input Parameters:
1238: + s - The PetscSection for the local field layout
1239: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1240: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1241: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1243: Output Parameter:
1244: . gsection - The PetscSection for the global field layout
1246: Note: This gives negative sizes and offsets to points not owned by this process
1248: Level: intermediate
1250: .seealso: PetscSectionCreate()
1251: @*/
1252: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1253: {
1254: PetscSection gs;
1255: const PetscInt *pind = NULL;
1256: PetscInt *recv = NULL, *neg = NULL;
1257: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1258: PetscErrorCode ierr;
1266: PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1267: PetscSectionGetChart(s, &pStart, &pEnd);
1268: PetscSectionSetChart(gs, pStart, pEnd);
1269: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1270: nlocal = nroots; /* The local/leaf space matches global/root space */
1271: /* Must allocate for all points visible to SF, which may be more than this section */
1272: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1273: PetscSFGetLeafRange(sf, NULL, &maxleaf);
1274: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1275: if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1276: PetscMalloc2(nroots,&neg,nlocal,&recv);
1277: PetscArrayzero(neg,nroots);
1278: }
1279: /* Mark all local points with negative dof */
1280: for (p = pStart; p < pEnd; ++p) {
1281: PetscSectionGetDof(s, p, &dof);
1282: PetscSectionSetDof(gs, p, dof);
1283: PetscSectionGetConstraintDof(s, p, &cdof);
1284: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1285: if (neg) neg[p] = -(dof+1);
1286: }
1287: PetscSectionSetUpBC(gs);
1288: if (gs->bcIndices) {PetscArraycpy(gs->bcIndices, s->bcIndices,gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]);}
1289: if (nroots >= 0) {
1290: PetscArrayzero(recv,nlocal);
1291: PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1292: PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1293: for (p = pStart; p < pEnd; ++p) {
1294: if (recv[p] < 0) {
1295: gs->atlasDof[p-pStart] = recv[p];
1296: PetscSectionGetDof(s, p, &dof);
1297: 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);
1298: }
1299: }
1300: }
1301: /* Calculate new sizes, get process offset, and calculate point offsets */
1302: if (s->perm) {ISGetIndices(s->perm, &pind);}
1303: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1304: const PetscInt q = pind ? pind[p] : p;
1306: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1307: gs->atlasOff[q] = off;
1308: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1309: }
1310: if (!localOffsets) {
1311: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1312: globalOff -= off;
1313: }
1314: for (p = pStart, off = 0; p < pEnd; ++p) {
1315: gs->atlasOff[p-pStart] += globalOff;
1316: if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1317: }
1318: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1319: /* Put in negative offsets for ghost points */
1320: if (nroots >= 0) {
1321: PetscArrayzero(recv,nlocal);
1322: PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1323: PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1324: for (p = pStart; p < pEnd; ++p) {
1325: if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1326: }
1327: }
1328: PetscFree2(neg,recv);
1329: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1330: *gsection = gs;
1331: return(0);
1332: }
1334: /*@
1335: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1336: the local section and an SF describing the section point overlap.
1338: Input Parameters:
1339: + s - The PetscSection for the local field layout
1340: . sf - The SF describing parallel layout of the section points
1341: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1342: . numExcludes - The number of exclusion ranges
1343: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1345: Output Parameter:
1346: . gsection - The PetscSection for the global field layout
1348: Note: This gives negative sizes and offsets to points not owned by this process
1350: Level: advanced
1352: .seealso: PetscSectionCreate()
1353: @*/
1354: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1355: {
1356: const PetscInt *pind = NULL;
1357: PetscInt *neg = NULL, *tmpOff = NULL;
1358: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1359: PetscErrorCode ierr;
1365: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1366: PetscSectionGetChart(s, &pStart, &pEnd);
1367: PetscSectionSetChart(*gsection, pStart, pEnd);
1368: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1369: if (nroots >= 0) {
1370: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1371: PetscCalloc1(nroots, &neg);
1372: if (nroots > pEnd-pStart) {
1373: PetscCalloc1(nroots, &tmpOff);
1374: } else {
1375: tmpOff = &(*gsection)->atlasDof[-pStart];
1376: }
1377: }
1378: /* Mark ghost points with negative dof */
1379: for (p = pStart; p < pEnd; ++p) {
1380: for (e = 0; e < numExcludes; ++e) {
1381: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1382: PetscSectionSetDof(*gsection, p, 0);
1383: break;
1384: }
1385: }
1386: if (e < numExcludes) continue;
1387: PetscSectionGetDof(s, p, &dof);
1388: PetscSectionSetDof(*gsection, p, dof);
1389: PetscSectionGetConstraintDof(s, p, &cdof);
1390: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1391: if (neg) neg[p] = -(dof+1);
1392: }
1393: PetscSectionSetUpBC(*gsection);
1394: if (nroots >= 0) {
1395: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1396: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1397: if (nroots > pEnd - pStart) {
1398: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1399: }
1400: }
1401: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1402: if (s->perm) {ISGetIndices(s->perm, &pind);}
1403: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1404: const PetscInt q = pind ? pind[p] : p;
1406: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1407: (*gsection)->atlasOff[q] = off;
1408: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1409: }
1410: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1411: globalOff -= off;
1412: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1413: (*gsection)->atlasOff[p] += globalOff;
1414: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1415: }
1416: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1417: /* Put in negative offsets for ghost points */
1418: if (nroots >= 0) {
1419: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1420: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1421: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1422: if (nroots > pEnd - pStart) {
1423: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1424: }
1425: }
1426: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1427: PetscFree(neg);
1428: return(0);
1429: }
1431: /*@
1432: PetscSectionGetPointLayout - Get the PetscLayout associated with the section points
1434: Collective on comm
1436: Input Parameters:
1437: + comm - The MPI_Comm
1438: - s - The PetscSection
1440: Output Parameter:
1441: . layout - The point layout for the section
1443: Note: This is usually called for the default global section.
1445: Level: advanced
1447: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1448: @*/
1449: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1450: {
1451: PetscInt pStart, pEnd, p, localSize = 0;
1455: PetscSectionGetChart(s, &pStart, &pEnd);
1456: for (p = pStart; p < pEnd; ++p) {
1457: PetscInt dof;
1459: PetscSectionGetDof(s, p, &dof);
1460: if (dof > 0) ++localSize;
1461: }
1462: PetscLayoutCreate(comm, layout);
1463: PetscLayoutSetLocalSize(*layout, localSize);
1464: PetscLayoutSetBlockSize(*layout, 1);
1465: PetscLayoutSetUp(*layout);
1466: return(0);
1467: }
1469: /*@
1470: PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.
1472: Collective on comm
1474: Input Parameters:
1475: + comm - The MPI_Comm
1476: - s - The PetscSection
1478: Output Parameter:
1479: . layout - The dof layout for the section
1481: Note: This is usually called for the default global section.
1483: Level: advanced
1485: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1486: @*/
1487: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1488: {
1489: PetscInt pStart, pEnd, p, localSize = 0;
1495: PetscSectionGetChart(s, &pStart, &pEnd);
1496: for (p = pStart; p < pEnd; ++p) {
1497: PetscInt dof,cdof;
1499: PetscSectionGetDof(s, p, &dof);
1500: PetscSectionGetConstraintDof(s, p, &cdof);
1501: if (dof-cdof > 0) localSize += dof-cdof;
1502: }
1503: PetscLayoutCreate(comm, layout);
1504: PetscLayoutSetLocalSize(*layout, localSize);
1505: PetscLayoutSetBlockSize(*layout, 1);
1506: PetscLayoutSetUp(*layout);
1507: return(0);
1508: }
1510: /*@
1511: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1513: Not collective
1515: Input Parameters:
1516: + s - the PetscSection
1517: - point - the point
1519: Output Parameter:
1520: . offset - the offset
1522: Level: intermediate
1524: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1525: @*/
1526: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1527: {
1531: if (PetscDefined(USE_DEBUG)) {
1532: 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);
1533: }
1534: *offset = s->atlasOff[point - s->pStart];
1535: return(0);
1536: }
1538: /*@
1539: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1541: Not collective
1543: Input Parameters:
1544: + s - the PetscSection
1545: . point - the point
1546: - offset - the offset
1548: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1550: Level: intermediate
1552: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1553: @*/
1554: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1555: {
1558: 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);
1559: s->atlasOff[point - s->pStart] = offset;
1560: return(0);
1561: }
1563: /*@
1564: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1566: Not collective
1568: Input Parameters:
1569: + s - the PetscSection
1570: . point - the point
1571: - field - the field
1573: Output Parameter:
1574: . offset - the offset
1576: Level: intermediate
1578: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1579: @*/
1580: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1581: {
1587: 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);
1588: PetscSectionGetOffset(s->field[field], point, offset);
1589: return(0);
1590: }
1592: /*@
1593: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given field at a point.
1595: Not collective
1597: Input Parameters:
1598: + s - the PetscSection
1599: . point - the point
1600: . field - the field
1601: - offset - the offset
1603: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1605: Level: intermediate
1607: .seealso: PetscSectionGetFieldOffset(), PetscSectionSetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1608: @*/
1609: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1610: {
1615: 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);
1616: PetscSectionSetOffset(s->field[field], point, offset);
1617: return(0);
1618: }
1620: /*@
1621: PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.
1623: Not collective
1625: Input Parameters:
1626: + s - the PetscSection
1627: . point - the point
1628: - field - the field
1630: Output Parameter:
1631: . offset - the offset
1633: Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1634: this point, what number is the first dof with this field.
1636: Level: advanced
1638: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1639: @*/
1640: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1641: {
1642: PetscInt off, foff;
1648: 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);
1649: PetscSectionGetOffset(s, point, &off);
1650: PetscSectionGetOffset(s->field[field], point, &foff);
1651: *offset = foff - off;
1652: return(0);
1653: }
1655: /*@
1656: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1658: Not collective
1660: Input Parameter:
1661: . s - the PetscSection
1663: Output Parameters:
1664: + start - the minimum offset
1665: - end - one more than the maximum offset
1667: Level: intermediate
1669: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1670: @*/
1671: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1672: {
1673: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1678: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1679: PetscSectionGetChart(s, &pStart, &pEnd);
1680: for (p = 0; p < pEnd-pStart; ++p) {
1681: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1683: if (off >= 0) {
1684: os = PetscMin(os, off);
1685: oe = PetscMax(oe, off+dof);
1686: }
1687: }
1688: if (start) *start = os;
1689: if (end) *end = oe;
1690: return(0);
1691: }
1693: /*@
1694: PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields
1696: Collective on s
1698: Input Parameter:
1699: + s - the PetscSection
1700: . len - the number of subfields
1701: - fields - the subfield numbers
1703: Output Parameter:
1704: . subs - the subsection
1706: Note: The section offsets now refer to a new, smaller vector.
1708: Level: advanced
1710: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1711: @*/
1712: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1713: {
1714: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
1718: if (!len) return(0);
1722: PetscSectionGetNumFields(s, &nF);
1723: if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1724: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1725: PetscSectionSetNumFields(*subs, len);
1726: for (f = 0; f < len; ++f) {
1727: const char *name = NULL;
1728: PetscInt numComp = 0;
1730: PetscSectionGetFieldName(s, fields[f], &name);
1731: PetscSectionSetFieldName(*subs, f, name);
1732: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1733: PetscSectionSetFieldComponents(*subs, f, numComp);
1734: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1735: PetscSectionGetComponentName(s, fields[f], c, &name);
1736: PetscSectionSetComponentName(*subs, f, c, name);
1737: }
1738: }
1739: PetscSectionGetChart(s, &pStart, &pEnd);
1740: PetscSectionSetChart(*subs, pStart, pEnd);
1741: for (p = pStart; p < pEnd; ++p) {
1742: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1744: for (f = 0; f < len; ++f) {
1745: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1746: PetscSectionSetFieldDof(*subs, p, f, fdof);
1747: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1748: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1749: dof += fdof;
1750: cdof += cfdof;
1751: }
1752: PetscSectionSetDof(*subs, p, dof);
1753: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1754: maxCdof = PetscMax(cdof, maxCdof);
1755: }
1756: PetscSectionSetUp(*subs);
1757: if (maxCdof) {
1758: PetscInt *indices;
1760: PetscMalloc1(maxCdof, &indices);
1761: for (p = pStart; p < pEnd; ++p) {
1762: PetscInt cdof;
1764: PetscSectionGetConstraintDof(*subs, p, &cdof);
1765: if (cdof) {
1766: const PetscInt *oldIndices = NULL;
1767: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1769: for (f = 0; f < len; ++f) {
1770: PetscInt oldFoff = 0, oldf;
1772: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1773: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1774: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1775: /* This can be sped up if we assume sorted fields */
1776: for (oldf = 0; oldf < fields[f]; ++oldf) {
1777: PetscInt oldfdof = 0;
1778: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1779: oldFoff += oldfdof;
1780: }
1781: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1782: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1783: numConst += cfdof;
1784: fOff += fdof;
1785: }
1786: PetscSectionSetConstraintIndices(*subs, p, indices);
1787: }
1788: }
1789: PetscFree(indices);
1790: }
1791: return(0);
1792: }
1794: /*@
1795: PetscSectionCreateSupersection - Create a new, larger section composed of the input sections
1797: Collective on s
1799: Input Parameters:
1800: + s - the input sections
1801: - len - the number of input sections
1803: Output Parameter:
1804: . supers - the supersection
1806: Note: The section offsets now refer to a new, larger vector.
1808: Level: advanced
1810: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1811: @*/
1812: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1813: {
1814: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1818: if (!len) return(0);
1819: for (i = 0; i < len; ++i) {
1820: PetscInt nf, pStarti, pEndi;
1822: PetscSectionGetNumFields(s[i], &nf);
1823: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1824: pStart = PetscMin(pStart, pStarti);
1825: pEnd = PetscMax(pEnd, pEndi);
1826: Nf += nf;
1827: }
1828: PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1829: PetscSectionSetNumFields(*supers, Nf);
1830: for (i = 0, f = 0; i < len; ++i) {
1831: PetscInt nf, fi, ci;
1833: PetscSectionGetNumFields(s[i], &nf);
1834: for (fi = 0; fi < nf; ++fi, ++f) {
1835: const char *name = NULL;
1836: PetscInt numComp = 0;
1838: PetscSectionGetFieldName(s[i], fi, &name);
1839: PetscSectionSetFieldName(*supers, f, name);
1840: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1841: PetscSectionSetFieldComponents(*supers, f, numComp);
1842: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1843: PetscSectionGetComponentName(s[i], fi, ci, &name);
1844: PetscSectionSetComponentName(*supers, f, ci, name);
1845: }
1846: }
1847: }
1848: PetscSectionSetChart(*supers, pStart, pEnd);
1849: for (p = pStart; p < pEnd; ++p) {
1850: PetscInt dof = 0, cdof = 0;
1852: for (i = 0, f = 0; i < len; ++i) {
1853: PetscInt nf, fi, pStarti, pEndi;
1854: PetscInt fdof = 0, cfdof = 0;
1856: PetscSectionGetNumFields(s[i], &nf);
1857: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1858: if ((p < pStarti) || (p >= pEndi)) continue;
1859: for (fi = 0; fi < nf; ++fi, ++f) {
1860: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1861: PetscSectionAddFieldDof(*supers, p, f, fdof);
1862: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1863: if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1864: dof += fdof;
1865: cdof += cfdof;
1866: }
1867: }
1868: PetscSectionSetDof(*supers, p, dof);
1869: if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1870: maxCdof = PetscMax(cdof, maxCdof);
1871: }
1872: PetscSectionSetUp(*supers);
1873: if (maxCdof) {
1874: PetscInt *indices;
1876: PetscMalloc1(maxCdof, &indices);
1877: for (p = pStart; p < pEnd; ++p) {
1878: PetscInt cdof;
1880: PetscSectionGetConstraintDof(*supers, p, &cdof);
1881: if (cdof) {
1882: PetscInt dof, numConst = 0, fOff = 0;
1884: for (i = 0, f = 0; i < len; ++i) {
1885: const PetscInt *oldIndices = NULL;
1886: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1888: PetscSectionGetNumFields(s[i], &nf);
1889: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1890: if ((p < pStarti) || (p >= pEndi)) continue;
1891: for (fi = 0; fi < nf; ++fi, ++f) {
1892: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1893: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1894: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1895: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1896: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1897: numConst += cfdof;
1898: }
1899: PetscSectionGetDof(s[i], p, &dof);
1900: fOff += dof;
1901: }
1902: PetscSectionSetConstraintIndices(*supers, p, indices);
1903: }
1904: }
1905: PetscFree(indices);
1906: }
1907: return(0);
1908: }
1910: /*@
1911: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
1913: Collective on s
1915: Input Parameters:
1916: + s - the PetscSection
1917: - subpointMap - a sorted list of points in the original mesh which are in the submesh
1919: Output Parameter:
1920: . subs - the subsection
1922: Note: The section offsets now refer to a new, smaller vector.
1924: Level: advanced
1926: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1927: @*/
1928: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1929: {
1930: const PetscInt *points = NULL, *indices = NULL;
1931: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;
1938: PetscSectionGetNumFields(s, &numFields);
1939: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1940: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1941: for (f = 0; f < numFields; ++f) {
1942: const char *name = NULL;
1943: PetscInt numComp = 0;
1945: PetscSectionGetFieldName(s, f, &name);
1946: PetscSectionSetFieldName(*subs, f, name);
1947: PetscSectionGetFieldComponents(s, f, &numComp);
1948: PetscSectionSetFieldComponents(*subs, f, numComp);
1949: for (c = 0; c < s->numFieldComponents[f]; ++c) {
1950: PetscSectionGetComponentName(s, f, c, &name);
1951: PetscSectionSetComponentName(*subs, f, c, name);
1952: }
1953: }
1954: /* For right now, we do not try to squeeze the subchart */
1955: if (subpointMap) {
1956: ISGetSize(subpointMap, &numSubpoints);
1957: ISGetIndices(subpointMap, &points);
1958: }
1959: PetscSectionGetChart(s, &pStart, &pEnd);
1960: PetscSectionSetChart(*subs, 0, numSubpoints);
1961: for (p = pStart; p < pEnd; ++p) {
1962: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1964: PetscFindInt(p, numSubpoints, points, &subp);
1965: if (subp < 0) continue;
1966: for (f = 0; f < numFields; ++f) {
1967: PetscSectionGetFieldDof(s, p, f, &fdof);
1968: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1969: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1970: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1971: }
1972: PetscSectionGetDof(s, p, &dof);
1973: PetscSectionSetDof(*subs, subp, dof);
1974: PetscSectionGetConstraintDof(s, p, &cdof);
1975: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1976: }
1977: PetscSectionSetUp(*subs);
1978: /* Change offsets to original offsets */
1979: for (p = pStart; p < pEnd; ++p) {
1980: PetscInt off, foff = 0;
1982: PetscFindInt(p, numSubpoints, points, &subp);
1983: if (subp < 0) continue;
1984: for (f = 0; f < numFields; ++f) {
1985: PetscSectionGetFieldOffset(s, p, f, &foff);
1986: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1987: }
1988: PetscSectionGetOffset(s, p, &off);
1989: PetscSectionSetOffset(*subs, subp, off);
1990: }
1991: /* Copy constraint indices */
1992: for (subp = 0; subp < numSubpoints; ++subp) {
1993: PetscInt cdof;
1995: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1996: if (cdof) {
1997: for (f = 0; f < numFields; ++f) {
1998: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1999: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2000: }
2001: PetscSectionGetConstraintIndices(s, points[subp], &indices);
2002: PetscSectionSetConstraintIndices(*subs, subp, indices);
2003: }
2004: }
2005: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
2006: return(0);
2007: }
2009: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2010: {
2011: PetscInt p;
2012: PetscMPIInt rank;
2016: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2017: PetscViewerASCIIPushSynchronized(viewer);
2018: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2019: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2020: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2021: PetscInt b;
2023: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2024: if (s->bcIndices) {
2025: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
2026: PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
2027: }
2028: }
2029: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2030: } else {
2031: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2032: }
2033: }
2034: PetscViewerFlush(viewer);
2035: PetscViewerASCIIPopSynchronized(viewer);
2036: if (s->sym) {
2037: PetscViewerASCIIPushTab(viewer);
2038: PetscSectionSymView(s->sym,viewer);
2039: PetscViewerASCIIPopTab(viewer);
2040: }
2041: return(0);
2042: }
2044: /*@C
2045: PetscSectionViewFromOptions - View from Options
2047: Collective on PetscSection
2049: Input Parameters:
2050: + A - the PetscSection object to view
2051: . obj - Optional object
2052: - name - command line option
2054: Level: intermediate
2055: .seealso: PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
2056: @*/
2057: PetscErrorCode PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
2058: {
2063: PetscObjectViewFromOptions((PetscObject)A,obj,name);
2064: return(0);
2065: }
2067: /*@C
2068: PetscSectionView - Views a PetscSection
2070: Collective on PetscSection
2072: Input Parameters:
2073: + s - the PetscSection object to view
2074: - v - the viewer
2076: Level: beginner
2078: .seealso PetscSectionCreate(), PetscSectionDestroy()
2079: @*/
2080: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2081: {
2082: PetscBool isascii;
2083: PetscInt f;
2088: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
2090: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
2091: if (isascii) {
2092: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
2093: if (s->numFields) {
2094: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
2095: for (f = 0; f < s->numFields; ++f) {
2096: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
2097: PetscSectionView_ASCII(s->field[f], viewer);
2098: }
2099: } else {
2100: PetscSectionView_ASCII(s, viewer);
2101: }
2102: }
2103: return(0);
2104: }
2106: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2107: {
2109: PetscSectionClosurePermVal clVal;
2112: if (!section->clHash) return(0);
2113: kh_foreach_value(section->clHash, clVal, {
2114: PetscFree(clVal.perm);
2115: PetscFree(clVal.invPerm);
2116: });
2117: kh_destroy(ClPerm, section->clHash);
2118: section->clHash = NULL;
2119: return(0);
2120: }
2122: /*@
2123: PetscSectionReset - Frees all section data.
2125: Not collective
2127: Input Parameters:
2128: . s - the PetscSection
2130: Level: beginner
2132: .seealso: PetscSection, PetscSectionCreate()
2133: @*/
2134: PetscErrorCode PetscSectionReset(PetscSection s)
2135: {
2136: PetscInt f, c;
2141: for (f = 0; f < s->numFields; ++f) {
2142: PetscSectionDestroy(&s->field[f]);
2143: PetscFree(s->fieldNames[f]);
2144: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2145: PetscFree(s->compNames[f][c]);
2146: }
2147: PetscFree(s->compNames[f]);
2148: }
2149: PetscFree(s->numFieldComponents);
2150: PetscFree(s->fieldNames);
2151: PetscFree(s->compNames);
2152: PetscFree(s->field);
2153: PetscSectionDestroy(&s->bc);
2154: PetscFree(s->bcIndices);
2155: PetscFree2(s->atlasDof, s->atlasOff);
2156: PetscSectionDestroy(&s->clSection);
2157: ISDestroy(&s->clPoints);
2158: ISDestroy(&s->perm);
2159: PetscSectionResetClosurePermutation(s);
2160: PetscSectionSymDestroy(&s->sym);
2161: PetscSectionDestroy(&s->clSection);
2162: ISDestroy(&s->clPoints);
2164: s->pStart = -1;
2165: s->pEnd = -1;
2166: s->maxDof = 0;
2167: s->setup = PETSC_FALSE;
2168: s->numFields = 0;
2169: s->clObj = NULL;
2170: return(0);
2171: }
2173: /*@
2174: PetscSectionDestroy - Frees a section object and frees its range if that exists.
2176: Not collective
2178: Input Parameters:
2179: . s - the PetscSection
2181: Level: beginner
2183: .seealso: PetscSection, PetscSectionCreate()
2184: @*/
2185: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2186: {
2190: if (!*s) return(0);
2192: if (--((PetscObject)(*s))->refct > 0) {
2193: *s = NULL;
2194: return(0);
2195: }
2196: PetscSectionReset(*s);
2197: PetscHeaderDestroy(s);
2198: return(0);
2199: }
2201: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2202: {
2203: const PetscInt p = point - s->pStart;
2207: *values = &baseArray[s->atlasOff[p]];
2208: return(0);
2209: }
2211: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2212: {
2213: PetscInt *array;
2214: const PetscInt p = point - s->pStart;
2215: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2216: PetscInt cDim = 0;
2221: PetscSectionGetConstraintDof(s, p, &cDim);
2222: array = &baseArray[s->atlasOff[p]];
2223: if (!cDim) {
2224: if (orientation >= 0) {
2225: const PetscInt dim = s->atlasDof[p];
2226: PetscInt i;
2228: if (mode == INSERT_VALUES) {
2229: for (i = 0; i < dim; ++i) array[i] = values[i];
2230: } else {
2231: for (i = 0; i < dim; ++i) array[i] += values[i];
2232: }
2233: } else {
2234: PetscInt offset = 0;
2235: PetscInt j = -1, field, i;
2237: for (field = 0; field < s->numFields; ++field) {
2238: const PetscInt dim = s->field[field]->atlasDof[p];
2240: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2241: offset += dim;
2242: }
2243: }
2244: } else {
2245: if (orientation >= 0) {
2246: const PetscInt dim = s->atlasDof[p];
2247: PetscInt cInd = 0, i;
2248: const PetscInt *cDof;
2250: PetscSectionGetConstraintIndices(s, point, &cDof);
2251: if (mode == INSERT_VALUES) {
2252: for (i = 0; i < dim; ++i) {
2253: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2254: array[i] = values[i];
2255: }
2256: } else {
2257: for (i = 0; i < dim; ++i) {
2258: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2259: array[i] += values[i];
2260: }
2261: }
2262: } else {
2263: const PetscInt *cDof;
2264: PetscInt offset = 0;
2265: PetscInt cOffset = 0;
2266: PetscInt j = 0, field;
2268: PetscSectionGetConstraintIndices(s, point, &cDof);
2269: for (field = 0; field < s->numFields; ++field) {
2270: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2271: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2272: const PetscInt sDim = dim - tDim;
2273: PetscInt cInd = 0, i,k;
2275: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2276: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2277: array[j] = values[k];
2278: }
2279: offset += dim;
2280: cOffset += dim - tDim;
2281: }
2282: }
2283: }
2284: return(0);
2285: }
2287: /*@C
2288: PetscSectionHasConstraints - Determine whether a section has constrained dofs
2290: Not collective
2292: Input Parameter:
2293: . s - The PetscSection
2295: Output Parameter:
2296: . hasConstraints - flag indicating that the section has constrained dofs
2298: Level: intermediate
2300: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2301: @*/
2302: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2303: {
2307: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2308: return(0);
2309: }
2311: /*@C
2312: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2314: Not collective
2316: Input Parameters:
2317: + s - The PetscSection
2318: - point - The point
2320: Output Parameter:
2321: . indices - The constrained dofs
2323: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
2325: Level: intermediate
2327: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2328: @*/
2329: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2330: {
2335: if (s->bc) {
2336: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2337: } else *indices = NULL;
2338: return(0);
2339: }
2341: /*@C
2342: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2344: Not collective
2346: Input Parameters:
2347: + s - The PetscSection
2348: . point - The point
2349: - indices - The constrained dofs
2351: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
2353: Level: intermediate
2355: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2356: @*/
2357: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2358: {
2363: if (s->bc) {
2364: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2365: }
2366: return(0);
2367: }
2369: /*@C
2370: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2372: Not collective
2374: Input Parameters:
2375: + s - The PetscSection
2376: . field - The field number
2377: - point - The point
2379: Output Parameter:
2380: . indices - The constrained dofs sorted in ascending order
2382: Notes:
2383: The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by PetscSectionGetConstraintDof().
2385: Fortran Note:
2386: In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()
2388: Level: intermediate
2390: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2391: @*/
2392: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2393: {
2398: 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);
2399: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2400: return(0);
2401: }
2403: /*@C
2404: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2406: Not collective
2408: Input Parameters:
2409: + s - The PetscSection
2410: . point - The point
2411: . field - The field number
2412: - indices - The constrained dofs
2414: Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()
2416: Level: intermediate
2418: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2419: @*/
2420: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2421: {
2426: 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);
2427: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2428: return(0);
2429: }
2431: /*@
2432: PetscSectionPermute - Reorder the section according to the input point permutation
2434: Collective on PetscSection
2436: Input Parameter:
2437: + section - The PetscSection object
2438: - perm - The point permutation, old point p becomes new point perm[p]
2440: Output Parameter:
2441: . sectionNew - The permuted PetscSection
2443: Level: intermediate
2445: .seealso: MatPermute()
2446: @*/
2447: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2448: {
2449: PetscSection s = section, sNew;
2450: const PetscInt *perm;
2451: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2452: PetscErrorCode ierr;
2458: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2459: PetscSectionGetNumFields(s, &numFields);
2460: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2461: for (f = 0; f < numFields; ++f) {
2462: const char *name;
2463: PetscInt numComp;
2465: PetscSectionGetFieldName(s, f, &name);
2466: PetscSectionSetFieldName(sNew, f, name);
2467: PetscSectionGetFieldComponents(s, f, &numComp);
2468: PetscSectionSetFieldComponents(sNew, f, numComp);
2469: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2470: PetscSectionGetComponentName(s, f, c, &name);
2471: PetscSectionSetComponentName(sNew, f, c, name);
2472: }
2473: }
2474: ISGetLocalSize(permutation, &numPoints);
2475: ISGetIndices(permutation, &perm);
2476: PetscSectionGetChart(s, &pStart, &pEnd);
2477: PetscSectionSetChart(sNew, pStart, pEnd);
2478: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2479: for (p = pStart; p < pEnd; ++p) {
2480: PetscInt dof, cdof;
2482: PetscSectionGetDof(s, p, &dof);
2483: PetscSectionSetDof(sNew, perm[p], dof);
2484: PetscSectionGetConstraintDof(s, p, &cdof);
2485: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2486: for (f = 0; f < numFields; ++f) {
2487: PetscSectionGetFieldDof(s, p, f, &dof);
2488: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2489: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2490: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2491: }
2492: }
2493: PetscSectionSetUp(sNew);
2494: for (p = pStart; p < pEnd; ++p) {
2495: const PetscInt *cind;
2496: PetscInt cdof;
2498: PetscSectionGetConstraintDof(s, p, &cdof);
2499: if (cdof) {
2500: PetscSectionGetConstraintIndices(s, p, &cind);
2501: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2502: }
2503: for (f = 0; f < numFields; ++f) {
2504: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2505: if (cdof) {
2506: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2507: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2508: }
2509: }
2510: }
2511: ISRestoreIndices(permutation, &perm);
2512: *sectionNew = sNew;
2513: return(0);
2514: }
2516: /*@
2517: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2519: Collective on section
2521: Input Parameters:
2522: + section - The PetscSection
2523: . obj - A PetscObject which serves as the key for this index
2524: . clSection - Section giving the size of the closure of each point
2525: - clPoints - IS giving the points in each closure
2527: Note: We compress out closure points with no dofs in this section
2529: Level: advanced
2531: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2532: @*/
2533: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2534: {
2541: if (section->clObj != obj) {PetscSectionResetClosurePermutation(section);}
2542: section->clObj = obj;
2543: PetscObjectReference((PetscObject)clSection);
2544: PetscObjectReference((PetscObject)clPoints);
2545: PetscSectionDestroy(§ion->clSection);
2546: ISDestroy(§ion->clPoints);
2547: section->clSection = clSection;
2548: section->clPoints = clPoints;
2549: return(0);
2550: }
2552: /*@
2553: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2555: Collective on section
2557: Input Parameters:
2558: + section - The PetscSection
2559: - obj - A PetscObject which serves as the key for this index
2561: Output Parameters:
2562: + clSection - Section giving the size of the closure of each point
2563: - clPoints - IS giving the points in each closure
2565: Note: We compress out closure points with no dofs in this section
2567: Level: advanced
2569: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2570: @*/
2571: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2572: {
2574: if (section->clObj == obj) {
2575: if (clSection) *clSection = section->clSection;
2576: if (clPoints) *clPoints = section->clPoints;
2577: } else {
2578: if (clSection) *clSection = NULL;
2579: if (clPoints) *clPoints = NULL;
2580: }
2581: return(0);
2582: }
2584: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2585: {
2586: PetscInt i;
2587: khiter_t iter;
2588: int new_entry;
2589: PetscSectionClosurePermKey key = {depth, clSize};
2590: PetscSectionClosurePermVal *val;
2594: if (section->clObj != obj) {
2595: PetscSectionDestroy(§ion->clSection);
2596: ISDestroy(§ion->clPoints);
2597: }
2598: section->clObj = obj;
2599: if (!section->clHash) {PetscClPermCreate(§ion->clHash);}
2600: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2601: val = &kh_val(section->clHash, iter);
2602: if (!new_entry) {
2603: PetscFree(val->perm);
2604: PetscFree(val->invPerm);
2605: }
2606: if (mode == PETSC_COPY_VALUES) {
2607: PetscMalloc1(clSize, &val->perm);
2608: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2609: PetscArraycpy(val->perm, clPerm, clSize);
2610: } else if (mode == PETSC_OWN_POINTER) {
2611: val->perm = clPerm;
2612: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2613: PetscMalloc1(clSize, &val->invPerm);
2614: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2615: return(0);
2616: }
2618: /*@
2619: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2621: Not Collective
2623: Input Parameters:
2624: + section - The PetscSection
2625: . obj - A PetscObject which serves as the key for this index (usually a DM)
2626: . depth - Depth of points on which to apply the given permutation
2627: - perm - Permutation of the cell dof closure
2629: Note:
2630: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2631: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2632: each topology and degree.
2634: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2635: exotic/enriched spaces on mixed topology meshes.
2637: Level: intermediate
2639: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2640: @*/
2641: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2642: {
2643: const PetscInt *clPerm = NULL;
2644: PetscInt clSize = 0;
2645: PetscErrorCode ierr;
2648: if (perm) {
2649: ISGetLocalSize(perm, &clSize);
2650: ISGetIndices(perm, &clPerm);
2651: }
2652: PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2653: if (perm) {ISRestoreIndices(perm, &clPerm);}
2654: return(0);
2655: }
2657: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2658: {
2662: if (section->clObj == obj) {
2663: PetscSectionClosurePermKey k = {depth, size};
2664: PetscSectionClosurePermVal v;
2665: PetscClPermGet(section->clHash, k, &v);
2666: if (perm) *perm = v.perm;
2667: } else {
2668: if (perm) *perm = NULL;
2669: }
2670: return(0);
2671: }
2673: /*@
2674: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2676: Not collective
2678: Input Parameters:
2679: + section - The PetscSection
2680: . obj - A PetscObject which serves as the key for this index (usually a DM)
2681: . depth - Depth stratum on which to obtain closure permutation
2682: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2684: Output Parameter:
2685: . perm - The dof closure permutation
2687: Note:
2688: The user must destroy the IS that is returned.
2690: Level: intermediate
2692: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2693: @*/
2694: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2695: {
2696: const PetscInt *clPerm;
2697: PetscErrorCode ierr;
2700: PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2701: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2702: return(0);
2703: }
2705: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2706: {
2710: if (section->clObj == obj && section->clHash) {
2711: PetscSectionClosurePermKey k = {depth, size};
2712: PetscSectionClosurePermVal v;
2713: PetscClPermGet(section->clHash, k, &v);
2714: if (perm) *perm = v.invPerm;
2715: } else {
2716: if (perm) *perm = NULL;
2717: }
2718: return(0);
2719: }
2721: /*@
2722: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2724: Not collective
2726: Input Parameters:
2727: + section - The PetscSection
2728: . obj - A PetscObject which serves as the key for this index (usually a DM)
2729: . depth - Depth stratum on which to obtain closure permutation
2730: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2732: Output Parameters:
2733: . perm - The dof closure permutation
2735: Note:
2736: The user must destroy the IS that is returned.
2738: Level: intermediate
2740: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2741: @*/
2742: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2743: {
2744: const PetscInt *clPerm;
2745: PetscErrorCode ierr;
2748: PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2749: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2750: return(0);
2751: }
2753: /*@
2754: PetscSectionGetField - Get the subsection associated with a single field
2756: Input Parameters:
2757: + s - The PetscSection
2758: - field - The field number
2760: Output Parameter:
2761: . subs - The subsection for the given field
2763: Level: intermediate
2765: .seealso: PetscSectionSetNumFields()
2766: @*/
2767: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2768: {
2772: 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);
2773: *subs = s->field[field];
2774: return(0);
2775: }
2777: PetscClassId PETSC_SECTION_SYM_CLASSID;
2778: PetscFunctionList PetscSectionSymList = NULL;
2780: /*@
2781: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
2783: Collective
2785: Input Parameter:
2786: . comm - the MPI communicator
2788: Output Parameter:
2789: . sym - pointer to the new set of symmetries
2791: Level: developer
2793: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2794: @*/
2795: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2796: {
2801: ISInitializePackage();
2802: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2803: return(0);
2804: }
2806: /*@C
2807: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
2809: Collective on PetscSectionSym
2811: Input Parameters:
2812: + sym - The section symmetry object
2813: - method - The name of the section symmetry type
2815: Level: developer
2817: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2818: @*/
2819: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2820: {
2821: PetscErrorCode (*r)(PetscSectionSym);
2822: PetscBool match;
2827: PetscObjectTypeCompare((PetscObject) sym, method, &match);
2828: if (match) return(0);
2830: PetscFunctionListFind(PetscSectionSymList,method,&r);
2831: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2832: if (sym->ops->destroy) {
2833: (*sym->ops->destroy)(sym);
2834: sym->ops->destroy = NULL;
2835: }
2836: (*r)(sym);
2837: PetscObjectChangeTypeName((PetscObject)sym,method);
2838: return(0);
2839: }
2842: /*@C
2843: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
2845: Not Collective
2847: Input Parameter:
2848: . sym - The section symmetry
2850: Output Parameter:
2851: . type - The index set type name
2853: Level: developer
2855: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2856: @*/
2857: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2858: {
2862: *type = ((PetscObject)sym)->type_name;
2863: return(0);
2864: }
2866: /*@C
2867: PetscSectionSymRegister - Adds a new section symmetry implementation
2869: Not Collective
2871: Input Parameters:
2872: + name - The name of a new user-defined creation routine
2873: - create_func - The creation routine itself
2875: Notes:
2876: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
2878: Level: developer
2880: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2881: @*/
2882: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2883: {
2887: ISInitializePackage();
2888: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2889: return(0);
2890: }
2892: /*@
2893: PetscSectionSymDestroy - Destroys a section symmetry.
2895: Collective on PetscSectionSym
2897: Input Parameters:
2898: . sym - the section symmetry
2900: Level: developer
2902: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2903: @*/
2904: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2905: {
2906: SymWorkLink link,next;
2910: if (!*sym) return(0);
2912: if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; return(0);}
2913: if ((*sym)->ops->destroy) {
2914: (*(*sym)->ops->destroy)(*sym);
2915: }
2916: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2917: for (link=(*sym)->workin; link; link=next) {
2918: next = link->next;
2919: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2920: PetscFree(link);
2921: }
2922: (*sym)->workin = NULL;
2923: PetscHeaderDestroy(sym);
2924: return(0);
2925: }
2927: /*@C
2928: PetscSectionSymView - Displays a section symmetry
2930: Collective on PetscSectionSym
2932: Input Parameters:
2933: + sym - the index set
2934: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
2936: Level: developer
2938: .seealso: PetscViewerASCIIOpen()
2939: @*/
2940: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2941: {
2946: if (!viewer) {
2947: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2948: }
2951: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2952: if (sym->ops->view) {
2953: (*sym->ops->view)(sym,viewer);
2954: }
2955: return(0);
2956: }
2958: /*@
2959: PetscSectionSetSym - Set the symmetries for the data referred to by the section
2961: Collective on PetscSection
2963: Input Parameters:
2964: + section - the section describing data layout
2965: - sym - the symmetry describing the affect of orientation on the access of the data
2967: Level: developer
2969: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2970: @*/
2971: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2972: {
2977: PetscSectionSymDestroy(&(section->sym));
2978: if (sym) {
2981: PetscObjectReference((PetscObject) sym);
2982: }
2983: section->sym = sym;
2984: return(0);
2985: }
2987: /*@
2988: PetscSectionGetSym - Get the symmetries for the data referred to by the section
2990: Not collective
2992: Input Parameters:
2993: . section - the section describing data layout
2995: Output Parameters:
2996: . sym - the symmetry describing the affect of orientation on the access of the data
2998: Level: developer
3000: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3001: @*/
3002: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3003: {
3006: *sym = section->sym;
3007: return(0);
3008: }
3010: /*@
3011: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3013: Collective on PetscSection
3015: Input Parameters:
3016: + section - the section describing data layout
3017: . field - the field number
3018: - sym - the symmetry describing the affect of orientation on the access of the data
3020: Level: developer
3022: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3023: @*/
3024: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3025: {
3030: 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);
3031: PetscSectionSetSym(section->field[field],sym);
3032: return(0);
3033: }
3035: /*@
3036: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3038: Collective on PetscSection
3040: Input Parameters:
3041: + section - the section describing data layout
3042: - field - the field number
3044: Output Parameters:
3045: . sym - the symmetry describing the affect of orientation on the access of the data
3047: Level: developer
3049: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3050: @*/
3051: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3052: {
3055: 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);
3056: *sym = section->field[field]->sym;
3057: return(0);
3058: }
3060: /*@C
3061: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
3063: Not collective
3065: Input Parameters:
3066: + section - the section
3067: . numPoints - the number of points
3068: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3069: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3070: context, see DMPlexGetConeOrientation()).
3072: Output Parameter:
3073: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3074: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3075: identity).
3077: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3078: .vb
3079: const PetscInt **perms;
3080: const PetscScalar **rots;
3081: PetscInt lOffset;
3083: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3084: for (i = 0, lOffset = 0; i < numPoints; i++) {
3085: PetscInt point = points[2*i], dof, sOffset;
3086: const PetscInt *perm = perms ? perms[i] : NULL;
3087: const PetscScalar *rot = rots ? rots[i] : NULL;
3089: PetscSectionGetDof(section,point,&dof);
3090: PetscSectionGetOffset(section,point,&sOffset);
3092: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3093: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3094: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3095: lOffset += dof;
3096: }
3097: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3098: .ve
3100: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3101: .vb
3102: const PetscInt **perms;
3103: const PetscScalar **rots;
3104: PetscInt lOffset;
3106: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3107: for (i = 0, lOffset = 0; i < numPoints; i++) {
3108: PetscInt point = points[2*i], dof, sOffset;
3109: const PetscInt *perm = perms ? perms[i] : NULL;
3110: const PetscScalar *rot = rots ? rots[i] : NULL;
3112: PetscSectionGetDof(section,point,&dof);
3113: PetscSectionGetOffset(section,point,&sOff);
3115: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3116: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3117: offset += dof;
3118: }
3119: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3120: .ve
3122: Level: developer
3124: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3125: @*/
3126: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3127: {
3128: PetscSectionSym sym;
3129: PetscErrorCode ierr;
3134: if (perms) *perms = NULL;
3135: if (rots) *rots = NULL;
3136: sym = section->sym;
3137: if (sym && (perms || rots)) {
3138: SymWorkLink link;
3140: if (sym->workin) {
3141: link = sym->workin;
3142: sym->workin = sym->workin->next;
3143: } else {
3144: PetscNewLog(sym,&link);
3145: }
3146: if (numPoints > link->numPoints) {
3147: PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3148: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3149: link->numPoints = numPoints;
3150: }
3151: link->next = sym->workout;
3152: sym->workout = link;
3153: PetscArrayzero((PetscInt**)link->perms,numPoints);
3154: PetscArrayzero((PetscInt**)link->rots,numPoints);
3155: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3156: if (perms) *perms = link->perms;
3157: if (rots) *rots = link->rots;
3158: }
3159: return(0);
3160: }
3162: /*@C
3163: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
3165: Not collective
3167: Input Parameters:
3168: + section - the section
3169: . numPoints - the number of points
3170: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3171: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3172: context, see DMPlexGetConeOrientation()).
3174: Output Parameter:
3175: + perms - The permutations for the given orientations: set to NULL at conclusion
3176: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3178: Level: developer
3180: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3181: @*/
3182: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3183: {
3184: PetscSectionSym sym;
3188: sym = section->sym;
3189: if (sym && (perms || rots)) {
3190: SymWorkLink *p,link;
3192: for (p=&sym->workout; (link=*p); p=&link->next) {
3193: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3194: *p = link->next;
3195: link->next = sym->workin;
3196: sym->workin = link;
3197: if (perms) *perms = NULL;
3198: if (rots) *rots = NULL;
3199: return(0);
3200: }
3201: }
3202: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3203: }
3204: return(0);
3205: }
3207: /*@C
3208: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
3210: Not collective
3212: Input Parameters:
3213: + section - the section
3214: . field - the field of the section
3215: . numPoints - the number of points
3216: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3217: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3218: context, see DMPlexGetConeOrientation()).
3220: Output Parameter:
3221: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3222: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3223: identity).
3225: Level: developer
3227: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3228: @*/
3229: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3230: {
3235: 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);
3236: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3237: return(0);
3238: }
3240: /*@C
3241: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
3243: Not collective
3245: Input Parameters:
3246: + section - the section
3247: . field - the field number
3248: . numPoints - the number of points
3249: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3250: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3251: context, see DMPlexGetConeOrientation()).
3253: Output Parameter:
3254: + perms - The permutations for the given orientations: set to NULL at conclusion
3255: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3257: Level: developer
3259: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3260: @*/
3261: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3262: {
3267: 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);
3268: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3269: return(0);
3270: }
3272: /*@
3273: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3275: Not collective
3277: Input Parameter:
3278: . s - the global PetscSection
3280: Output Parameters:
3281: . flg - the flag
3283: Level: developer
3285: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3286: @*/
3287: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3288: {
3291: *flg = s->useFieldOff;
3292: return(0);
3293: }
3295: /*@
3296: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3298: Not collective
3300: Input Parameters:
3301: + s - the global PetscSection
3302: - flg - the flag
3304: Level: developer
3306: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3307: @*/
3308: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3309: {
3312: s->useFieldOff = flg;
3313: return(0);
3314: }
3316: #define PetscSectionExpandPoints_Loop(TYPE) \
3317: { \
3318: PetscInt i, n, o0, o1, size; \
3319: TYPE *a0 = (TYPE*)origArray, *a1; \
3320: PetscSectionGetStorageSize(s, &size); \
3321: PetscMalloc1(size, &a1); \
3322: for (i=0; i<npoints; i++) { \
3323: PetscSectionGetOffset(origSection, points_[i], &o0); \
3324: PetscSectionGetOffset(s, i, &o1); \
3325: PetscSectionGetDof(s, i, &n); \
3326: PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3327: } \
3328: *newArray = (void*)a1; \
3329: }
3331: /*@
3332: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3334: Not collective
3336: Input Parameters:
3337: + origSection - the PetscSection describing the layout of the array
3338: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3339: . origArray - the array; its size must be equal to the storage size of origSection
3340: - points - IS with points to extract; its indices must lie in the chart of origSection
3342: Output Parameters:
3343: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3344: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3346: Level: developer
3348: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3349: @*/
3350: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3351: {
3352: PetscSection s;
3353: const PetscInt *points_;
3354: PetscInt i, n, npoints, pStart, pEnd;
3355: PetscMPIInt unitsize;
3356: PetscErrorCode ierr;
3364: MPI_Type_size(dataType, &unitsize);
3365: ISGetLocalSize(points, &npoints);
3366: ISGetIndices(points, &points_);
3367: PetscSectionGetChart(origSection, &pStart, &pEnd);
3368: PetscSectionCreate(PETSC_COMM_SELF, &s);
3369: PetscSectionSetChart(s, 0, npoints);
3370: for (i=0; i<npoints; i++) {
3371: if (PetscUnlikely(points_[i] < pStart || points_[i] >= pEnd)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %d (index %d) in input IS out of input section's chart", points_[i], i);
3372: PetscSectionGetDof(origSection, points_[i], &n);
3373: PetscSectionSetDof(s, i, n);
3374: }
3375: PetscSectionSetUp(s);
3376: if (newArray) {
3377: if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);}
3378: else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);}
3379: else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);}
3380: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3381: }
3382: if (newSection) {
3383: *newSection = s;
3384: } else {
3385: PetscSectionDestroy(&s);
3386: }
3387: return(0);
3388: }