Actual source code: section.c
petsc-3.14.6 2021-03-30
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);
1292: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
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);
1323: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
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);
1396: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
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);
1421: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
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 caleld 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 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: PetscSectionGetOffset(), 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: PetscFree(s->compNames[f]);
2147: }
2148: PetscFree(s->numFieldComponents);
2149: PetscFree(s->fieldNames);
2150: PetscFree(s->compNames);
2151: PetscFree(s->field);
2152: PetscSectionDestroy(&s->bc);
2153: PetscFree(s->bcIndices);
2154: PetscFree2(s->atlasDof, s->atlasOff);
2155: PetscSectionDestroy(&s->clSection);
2156: ISDestroy(&s->clPoints);
2157: ISDestroy(&s->perm);
2158: PetscSectionResetClosurePermutation(s);
2159: PetscSectionSymDestroy(&s->sym);
2160: PetscSectionDestroy(&s->clSection);
2161: ISDestroy(&s->clPoints);
2163: s->pStart = -1;
2164: s->pEnd = -1;
2165: s->maxDof = 0;
2166: s->setup = PETSC_FALSE;
2167: s->numFields = 0;
2168: s->clObj = NULL;
2169: return(0);
2170: }
2172: /*@
2173: PetscSectionDestroy - Frees a section object and frees its range if that exists.
2175: Not collective
2177: Input Parameters:
2178: . s - the PetscSection
2180: Level: beginner
2182: .seealso: PetscSection, PetscSectionCreate()
2183: @*/
2184: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2185: {
2189: if (!*s) return(0);
2191: if (--((PetscObject)(*s))->refct > 0) {
2192: *s = NULL;
2193: return(0);
2194: }
2195: PetscSectionReset(*s);
2196: PetscHeaderDestroy(s);
2197: return(0);
2198: }
2200: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2201: {
2202: const PetscInt p = point - s->pStart;
2206: *values = &baseArray[s->atlasOff[p]];
2207: return(0);
2208: }
2210: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2211: {
2212: PetscInt *array;
2213: const PetscInt p = point - s->pStart;
2214: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2215: PetscInt cDim = 0;
2220: PetscSectionGetConstraintDof(s, p, &cDim);
2221: array = &baseArray[s->atlasOff[p]];
2222: if (!cDim) {
2223: if (orientation >= 0) {
2224: const PetscInt dim = s->atlasDof[p];
2225: PetscInt i;
2227: if (mode == INSERT_VALUES) {
2228: for (i = 0; i < dim; ++i) array[i] = values[i];
2229: } else {
2230: for (i = 0; i < dim; ++i) array[i] += values[i];
2231: }
2232: } else {
2233: PetscInt offset = 0;
2234: PetscInt j = -1, field, i;
2236: for (field = 0; field < s->numFields; ++field) {
2237: const PetscInt dim = s->field[field]->atlasDof[p];
2239: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2240: offset += dim;
2241: }
2242: }
2243: } else {
2244: if (orientation >= 0) {
2245: const PetscInt dim = s->atlasDof[p];
2246: PetscInt cInd = 0, i;
2247: const PetscInt *cDof;
2249: PetscSectionGetConstraintIndices(s, point, &cDof);
2250: if (mode == INSERT_VALUES) {
2251: for (i = 0; i < dim; ++i) {
2252: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2253: array[i] = values[i];
2254: }
2255: } else {
2256: for (i = 0; i < dim; ++i) {
2257: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2258: array[i] += values[i];
2259: }
2260: }
2261: } else {
2262: const PetscInt *cDof;
2263: PetscInt offset = 0;
2264: PetscInt cOffset = 0;
2265: PetscInt j = 0, field;
2267: PetscSectionGetConstraintIndices(s, point, &cDof);
2268: for (field = 0; field < s->numFields; ++field) {
2269: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2270: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2271: const PetscInt sDim = dim - tDim;
2272: PetscInt cInd = 0, i,k;
2274: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2275: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2276: array[j] = values[k];
2277: }
2278: offset += dim;
2279: cOffset += dim - tDim;
2280: }
2281: }
2282: }
2283: return(0);
2284: }
2286: /*@C
2287: PetscSectionHasConstraints - Determine whether a section has constrained dofs
2289: Not collective
2291: Input Parameter:
2292: . s - The PetscSection
2294: Output Parameter:
2295: . hasConstraints - flag indicating that the section has constrained dofs
2297: Level: intermediate
2299: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2300: @*/
2301: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2302: {
2306: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2307: return(0);
2308: }
2310: /*@C
2311: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2313: Not collective
2315: Input Parameters:
2316: + s - The PetscSection
2317: - point - The point
2319: Output Parameter:
2320: . indices - The constrained dofs
2322: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
2324: Level: intermediate
2326: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2327: @*/
2328: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2329: {
2334: if (s->bc) {
2335: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2336: } else *indices = NULL;
2337: return(0);
2338: }
2340: /*@C
2341: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2343: Not collective
2345: Input Parameters:
2346: + s - The PetscSection
2347: . point - The point
2348: - indices - The constrained dofs
2350: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
2352: Level: intermediate
2354: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2355: @*/
2356: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2357: {
2362: if (s->bc) {
2363: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2364: }
2365: return(0);
2366: }
2368: /*@C
2369: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2371: Not collective
2373: Input Parameters:
2374: + s - The PetscSection
2375: . field - The field number
2376: - point - The point
2378: Output Parameter:
2379: . indices - The constrained dofs sorted in ascending order
2381: Notes:
2382: 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().
2384: Fortran Note:
2385: In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()
2387: Level: intermediate
2389: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2390: @*/
2391: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2392: {
2397: 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);
2398: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2399: return(0);
2400: }
2402: /*@C
2403: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2405: Not collective
2407: Input Parameters:
2408: + s - The PetscSection
2409: . point - The point
2410: . field - The field number
2411: - indices - The constrained dofs
2413: Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()
2415: Level: intermediate
2417: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2418: @*/
2419: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2420: {
2425: 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);
2426: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2427: return(0);
2428: }
2430: /*@
2431: PetscSectionPermute - Reorder the section according to the input point permutation
2433: Collective on PetscSection
2435: Input Parameter:
2436: + section - The PetscSection object
2437: - perm - The point permutation, old point p becomes new point perm[p]
2439: Output Parameter:
2440: . sectionNew - The permuted PetscSection
2442: Level: intermediate
2444: .seealso: MatPermute()
2445: @*/
2446: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2447: {
2448: PetscSection s = section, sNew;
2449: const PetscInt *perm;
2450: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2451: PetscErrorCode ierr;
2457: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2458: PetscSectionGetNumFields(s, &numFields);
2459: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2460: for (f = 0; f < numFields; ++f) {
2461: const char *name;
2462: PetscInt numComp;
2464: PetscSectionGetFieldName(s, f, &name);
2465: PetscSectionSetFieldName(sNew, f, name);
2466: PetscSectionGetFieldComponents(s, f, &numComp);
2467: PetscSectionSetFieldComponents(sNew, f, numComp);
2468: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2469: PetscSectionGetComponentName(s, f, c, &name);
2470: PetscSectionSetComponentName(sNew, f, c, name);
2471: }
2472: }
2473: ISGetLocalSize(permutation, &numPoints);
2474: ISGetIndices(permutation, &perm);
2475: PetscSectionGetChart(s, &pStart, &pEnd);
2476: PetscSectionSetChart(sNew, pStart, pEnd);
2477: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2478: for (p = pStart; p < pEnd; ++p) {
2479: PetscInt dof, cdof;
2481: PetscSectionGetDof(s, p, &dof);
2482: PetscSectionSetDof(sNew, perm[p], dof);
2483: PetscSectionGetConstraintDof(s, p, &cdof);
2484: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2485: for (f = 0; f < numFields; ++f) {
2486: PetscSectionGetFieldDof(s, p, f, &dof);
2487: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2488: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2489: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2490: }
2491: }
2492: PetscSectionSetUp(sNew);
2493: for (p = pStart; p < pEnd; ++p) {
2494: const PetscInt *cind;
2495: PetscInt cdof;
2497: PetscSectionGetConstraintDof(s, p, &cdof);
2498: if (cdof) {
2499: PetscSectionGetConstraintIndices(s, p, &cind);
2500: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2501: }
2502: for (f = 0; f < numFields; ++f) {
2503: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2504: if (cdof) {
2505: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2506: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2507: }
2508: }
2509: }
2510: ISRestoreIndices(permutation, &perm);
2511: *sectionNew = sNew;
2512: return(0);
2513: }
2515: /* TODO: the next three functions should be moved to sf/utils */
2516: #include <petsc/private/sfimpl.h>
2518: /*@C
2519: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
2521: Collective on sf
2523: Input Parameters:
2524: + sf - The SF
2525: - rootSection - Section defined on root space
2527: Output Parameters:
2528: + remoteOffsets - root offsets in leaf storage, or NULL
2529: - leafSection - Section defined on the leaf space
2531: Level: advanced
2533: .seealso: PetscSFCreate()
2534: @*/
2535: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2536: {
2537: PetscSF embedSF;
2538: const PetscInt *indices;
2539: IS selected;
2540: PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
2541: PetscBool *sub, hasc;
2545: PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
2546: PetscSectionGetNumFields(rootSection, &numFields);
2547: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2548: PetscMalloc1(numFields+2, &sub);
2549: sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
2550: for (f = 0; f < numFields; ++f) {
2551: PetscSectionSym sym;
2552: const char *name = NULL;
2553: PetscInt numComp = 0;
2555: sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
2556: PetscSectionGetFieldComponents(rootSection, f, &numComp);
2557: PetscSectionGetFieldName(rootSection, f, &name);
2558: PetscSectionGetFieldSym(rootSection, f, &sym);
2559: PetscSectionSetFieldComponents(leafSection, f, numComp);
2560: PetscSectionSetFieldName(leafSection, f, name);
2561: PetscSectionSetFieldSym(leafSection, f, sym);
2562: for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
2563: PetscSectionGetComponentName(rootSection, f, c, &name);
2564: PetscSectionSetComponentName(leafSection, f, c, name);
2565: }
2566: }
2567: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2568: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2569: rpEnd = PetscMin(rpEnd,nroots);
2570: rpEnd = PetscMax(rpStart,rpEnd);
2571: /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
2572: sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
2573: MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
2574: if (sub[0]) {
2575: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2576: ISGetIndices(selected, &indices);
2577: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2578: ISRestoreIndices(selected, &indices);
2579: ISDestroy(&selected);
2580: } else {
2581: PetscObjectReference((PetscObject)sf);
2582: embedSF = sf;
2583: }
2584: PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
2585: lpEnd++;
2587: PetscSectionSetChart(leafSection, lpStart, lpEnd);
2589: /* Constrained dof section */
2590: hasc = sub[1];
2591: for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);
2593: /* Could fuse these at the cost of copies and extra allocation */
2594: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2595: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2596: if (sub[1]) {
2597: PetscSectionCheckConstraints_Static(rootSection);
2598: PetscSectionCheckConstraints_Static(leafSection);
2599: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2600: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2601: }
2602: for (f = 0; f < numFields; ++f) {
2603: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2604: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2605: if (sub[2+f]) {
2606: PetscSectionCheckConstraints_Static(rootSection->field[f]);
2607: PetscSectionCheckConstraints_Static(leafSection->field[f]);
2608: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2609: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2610: }
2611: }
2612: if (remoteOffsets) {
2613: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2614: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2615: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2616: }
2617: PetscSectionSetUp(leafSection);
2618: if (hasc) { /* need to communicate bcIndices */
2619: PetscSF bcSF;
2620: PetscInt *rOffBc;
2622: PetscMalloc1(lpEnd - lpStart, &rOffBc);
2623: if (sub[1]) {
2624: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2625: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2626: PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
2627: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2628: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2629: PetscSFDestroy(&bcSF);
2630: }
2631: for (f = 0; f < numFields; ++f) {
2632: if (sub[2+f]) {
2633: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2634: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2635: PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
2636: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2637: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2638: PetscSFDestroy(&bcSF);
2639: }
2640: }
2641: PetscFree(rOffBc);
2642: }
2643: PetscSFDestroy(&embedSF);
2644: PetscFree(sub);
2645: PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
2646: return(0);
2647: }
2649: /*@C
2650: PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
2652: Collective on sf
2654: Input Parameters:
2655: + sf - The SF
2656: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
2657: - leafSection - Data layout of local points for incoming data (this is layout for SF leaves)
2659: Output Parameter:
2660: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2662: Level: developer
2664: .seealso: PetscSFCreate()
2665: @*/
2666: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2667: {
2668: PetscSF embedSF;
2669: const PetscInt *indices;
2670: IS selected;
2671: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2672: PetscErrorCode ierr;
2675: *remoteOffsets = NULL;
2676: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2677: if (numRoots < 0) return(0);
2678: PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
2679: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2680: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2681: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2682: ISGetIndices(selected, &indices);
2683: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2684: ISRestoreIndices(selected, &indices);
2685: ISDestroy(&selected);
2686: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2687: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2688: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2689: PetscSFDestroy(&embedSF);
2690: PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
2691: return(0);
2692: }
2694: /*@C
2695: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2697: Collective on sf
2699: Input Parameters:
2700: + sf - The SF
2701: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2702: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2703: - leafSection - Data layout of local points for incoming data (this is the distributed section)
2705: Output Parameters:
2706: - sectionSF - The new SF
2708: Note: Either rootSection or remoteOffsets can be specified
2710: Level: advanced
2712: .seealso: PetscSFCreate()
2713: @*/
2714: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2715: {
2716: MPI_Comm comm;
2717: const PetscInt *localPoints;
2718: const PetscSFNode *remotePoints;
2719: PetscInt lpStart, lpEnd;
2720: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2721: PetscInt *localIndices;
2722: PetscSFNode *remoteIndices;
2723: PetscInt i, ind;
2724: PetscErrorCode ierr;
2732: PetscObjectGetComm((PetscObject)sf,&comm);
2733: PetscSFCreate(comm, sectionSF);
2734: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2735: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2736: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2737: if (numRoots < 0) return(0);
2738: PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
2739: for (i = 0; i < numPoints; ++i) {
2740: PetscInt localPoint = localPoints ? localPoints[i] : i;
2741: PetscInt dof;
2743: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2744: PetscSectionGetDof(leafSection, localPoint, &dof);
2745: numIndices += dof;
2746: }
2747: }
2748: PetscMalloc1(numIndices, &localIndices);
2749: PetscMalloc1(numIndices, &remoteIndices);
2750: /* Create new index graph */
2751: for (i = 0, ind = 0; i < numPoints; ++i) {
2752: PetscInt localPoint = localPoints ? localPoints[i] : i;
2753: PetscInt rank = remotePoints[i].rank;
2755: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2756: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2757: PetscInt loff, dof, d;
2759: PetscSectionGetOffset(leafSection, localPoint, &loff);
2760: PetscSectionGetDof(leafSection, localPoint, &dof);
2761: for (d = 0; d < dof; ++d, ++ind) {
2762: localIndices[ind] = loff+d;
2763: remoteIndices[ind].rank = rank;
2764: remoteIndices[ind].index = remoteOffset+d;
2765: }
2766: }
2767: }
2768: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2769: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2770: PetscSFSetUp(*sectionSF);
2771: PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
2772: return(0);
2773: }
2775: /*@
2776: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2778: Collective on section
2780: Input Parameters:
2781: + section - The PetscSection
2782: . obj - A PetscObject which serves as the key for this index
2783: . clSection - Section giving the size of the closure of each point
2784: - clPoints - IS giving the points in each closure
2786: Note: We compress out closure points with no dofs in this section
2788: Level: advanced
2790: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2791: @*/
2792: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2793: {
2800: if (section->clObj != obj) {PetscSectionResetClosurePermutation(section);}
2801: section->clObj = obj;
2802: PetscObjectReference((PetscObject)clSection);
2803: PetscObjectReference((PetscObject)clPoints);
2804: PetscSectionDestroy(§ion->clSection);
2805: ISDestroy(§ion->clPoints);
2806: section->clSection = clSection;
2807: section->clPoints = clPoints;
2808: return(0);
2809: }
2811: /*@
2812: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2814: Collective on section
2816: Input Parameters:
2817: + section - The PetscSection
2818: - obj - A PetscObject which serves as the key for this index
2820: Output Parameters:
2821: + clSection - Section giving the size of the closure of each point
2822: - clPoints - IS giving the points in each closure
2824: Note: We compress out closure points with no dofs in this section
2826: Level: advanced
2828: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2829: @*/
2830: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2831: {
2833: if (section->clObj == obj) {
2834: if (clSection) *clSection = section->clSection;
2835: if (clPoints) *clPoints = section->clPoints;
2836: } else {
2837: if (clSection) *clSection = NULL;
2838: if (clPoints) *clPoints = NULL;
2839: }
2840: return(0);
2841: }
2843: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2844: {
2845: PetscInt i;
2846: khiter_t iter;
2847: int new_entry;
2848: PetscSectionClosurePermKey key = {depth, clSize};
2849: PetscSectionClosurePermVal *val;
2853: if (section->clObj != obj) {
2854: PetscSectionDestroy(§ion->clSection);
2855: ISDestroy(§ion->clPoints);
2856: }
2857: section->clObj = obj;
2858: if (!section->clHash) {PetscClPermCreate(§ion->clHash);}
2859: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2860: val = &kh_val(section->clHash, iter);
2861: if (!new_entry) {
2862: PetscFree(val->perm);
2863: PetscFree(val->invPerm);
2864: }
2865: if (mode == PETSC_COPY_VALUES) {
2866: PetscMalloc1(clSize, &val->perm);
2867: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2868: PetscArraycpy(val->perm, clPerm, clSize);
2869: } else if (mode == PETSC_OWN_POINTER) {
2870: val->perm = clPerm;
2871: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2872: PetscMalloc1(clSize, &val->invPerm);
2873: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2874: return(0);
2875: }
2877: /*@
2878: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2880: Not Collective
2882: Input Parameters:
2883: + section - The PetscSection
2884: . obj - A PetscObject which serves as the key for this index (usually a DM)
2885: . depth - Depth of points on which to apply the given permutation
2886: - perm - Permutation of the cell dof closure
2888: Note:
2889: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2890: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2891: each topology and degree.
2893: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2894: exotic/enriched spaces on mixed topology meshes.
2896: Level: intermediate
2898: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2899: @*/
2900: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2901: {
2902: const PetscInt *clPerm = NULL;
2903: PetscInt clSize = 0;
2904: PetscErrorCode ierr;
2907: if (perm) {
2908: ISGetLocalSize(perm, &clSize);
2909: ISGetIndices(perm, &clPerm);
2910: }
2911: PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2912: if (perm) {ISRestoreIndices(perm, &clPerm);}
2913: return(0);
2914: }
2916: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2917: {
2921: if (section->clObj == obj) {
2922: PetscSectionClosurePermKey k = {depth, size};
2923: PetscSectionClosurePermVal v;
2924: PetscClPermGet(section->clHash, k, &v);
2925: if (perm) *perm = v.perm;
2926: } else {
2927: if (perm) *perm = NULL;
2928: }
2929: return(0);
2930: }
2932: /*@
2933: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2935: Not collective
2937: Input Parameters:
2938: + section - The PetscSection
2939: . obj - A PetscObject which serves as the key for this index (usually a DM)
2940: . depth - Depth stratum on which to obtain closure permutation
2941: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2943: Output Parameter:
2944: . perm - The dof closure permutation
2946: Note:
2947: The user must destroy the IS that is returned.
2949: Level: intermediate
2951: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2952: @*/
2953: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2954: {
2955: const PetscInt *clPerm;
2956: PetscErrorCode ierr;
2959: PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2960: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2961: return(0);
2962: }
2964: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2965: {
2969: if (section->clObj == obj && section->clHash) {
2970: PetscSectionClosurePermKey k = {depth, size};
2971: PetscSectionClosurePermVal v;
2972: PetscClPermGet(section->clHash, k, &v);
2973: if (perm) *perm = v.invPerm;
2974: } else {
2975: if (perm) *perm = NULL;
2976: }
2977: return(0);
2978: }
2980: /*@
2981: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2983: Not collective
2985: Input Parameters:
2986: + section - The PetscSection
2987: . obj - A PetscObject which serves as the key for this index (usually a DM)
2988: . depth - Depth stratum on which to obtain closure permutation
2989: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2991: Output Parameters:
2992: . perm - The dof closure permutation
2994: Note:
2995: The user must destroy the IS that is returned.
2997: Level: intermediate
2999: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
3000: @*/
3001: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3002: {
3003: const PetscInt *clPerm;
3004: PetscErrorCode ierr;
3007: PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
3008: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
3009: return(0);
3010: }
3012: /*@
3013: PetscSectionGetField - Get the subsection associated with a single field
3015: Input Parameters:
3016: + s - The PetscSection
3017: - field - The field number
3019: Output Parameter:
3020: . subs - The subsection for the given field
3022: Level: intermediate
3024: .seealso: PetscSectionSetNumFields()
3025: @*/
3026: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3027: {
3031: 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);
3032: *subs = s->field[field];
3033: return(0);
3034: }
3036: PetscClassId PETSC_SECTION_SYM_CLASSID;
3037: PetscFunctionList PetscSectionSymList = NULL;
3039: /*@
3040: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
3042: Collective
3044: Input Parameter:
3045: . comm - the MPI communicator
3047: Output Parameter:
3048: . sym - pointer to the new set of symmetries
3050: Level: developer
3052: .seealso: PetscSectionSym, PetscSectionSymDestroy()
3053: @*/
3054: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3055: {
3060: ISInitializePackage();
3061: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
3062: return(0);
3063: }
3065: /*@C
3066: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
3068: Collective on PetscSectionSym
3070: Input Parameters:
3071: + sym - The section symmetry object
3072: - method - The name of the section symmetry type
3074: Level: developer
3076: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
3077: @*/
3078: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3079: {
3080: PetscErrorCode (*r)(PetscSectionSym);
3081: PetscBool match;
3086: PetscObjectTypeCompare((PetscObject) sym, method, &match);
3087: if (match) return(0);
3089: PetscFunctionListFind(PetscSectionSymList,method,&r);
3090: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3091: if (sym->ops->destroy) {
3092: (*sym->ops->destroy)(sym);
3093: sym->ops->destroy = NULL;
3094: }
3095: (*r)(sym);
3096: PetscObjectChangeTypeName((PetscObject)sym,method);
3097: return(0);
3098: }
3101: /*@C
3102: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
3104: Not Collective
3106: Input Parameter:
3107: . sym - The section symmetry
3109: Output Parameter:
3110: . type - The index set type name
3112: Level: developer
3114: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
3115: @*/
3116: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3117: {
3121: *type = ((PetscObject)sym)->type_name;
3122: return(0);
3123: }
3125: /*@C
3126: PetscSectionSymRegister - Adds a new section symmetry implementation
3128: Not Collective
3130: Input Parameters:
3131: + name - The name of a new user-defined creation routine
3132: - create_func - The creation routine itself
3134: Notes:
3135: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
3137: Level: developer
3139: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
3140: @*/
3141: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3142: {
3146: ISInitializePackage();
3147: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
3148: return(0);
3149: }
3151: /*@
3152: PetscSectionSymDestroy - Destroys a section symmetry.
3154: Collective on PetscSectionSym
3156: Input Parameters:
3157: . sym - the section symmetry
3159: Level: developer
3161: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
3162: @*/
3163: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3164: {
3165: SymWorkLink link,next;
3169: if (!*sym) return(0);
3171: if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; return(0);}
3172: if ((*sym)->ops->destroy) {
3173: (*(*sym)->ops->destroy)(*sym);
3174: }
3175: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3176: for (link=(*sym)->workin; link; link=next) {
3177: next = link->next;
3178: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
3179: PetscFree(link);
3180: }
3181: (*sym)->workin = NULL;
3182: PetscHeaderDestroy(sym);
3183: return(0);
3184: }
3186: /*@C
3187: PetscSectionSymView - Displays a section symmetry
3189: Collective on PetscSectionSym
3191: Input Parameters:
3192: + sym - the index set
3193: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
3195: Level: developer
3197: .seealso: PetscViewerASCIIOpen()
3198: @*/
3199: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3200: {
3205: if (!viewer) {
3206: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
3207: }
3210: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3211: if (sym->ops->view) {
3212: (*sym->ops->view)(sym,viewer);
3213: }
3214: return(0);
3215: }
3217: /*@
3218: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3220: Collective on PetscSection
3222: Input Parameters:
3223: + section - the section describing data layout
3224: - sym - the symmetry describing the affect of orientation on the access of the data
3226: Level: developer
3228: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3229: @*/
3230: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3231: {
3236: PetscSectionSymDestroy(&(section->sym));
3237: if (sym) {
3240: PetscObjectReference((PetscObject) sym);
3241: }
3242: section->sym = sym;
3243: return(0);
3244: }
3246: /*@
3247: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3249: Not collective
3251: Input Parameters:
3252: . section - the section describing data layout
3254: Output Parameters:
3255: . sym - the symmetry describing the affect of orientation on the access of the data
3257: Level: developer
3259: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3260: @*/
3261: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3262: {
3265: *sym = section->sym;
3266: return(0);
3267: }
3269: /*@
3270: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3272: Collective on PetscSection
3274: Input Parameters:
3275: + section - the section describing data layout
3276: . field - the field number
3277: - sym - the symmetry describing the affect of orientation on the access of the data
3279: Level: developer
3281: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3282: @*/
3283: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3284: {
3289: 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);
3290: PetscSectionSetSym(section->field[field],sym);
3291: return(0);
3292: }
3294: /*@
3295: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3297: Collective on PetscSection
3299: Input Parameters:
3300: + section - the section describing data layout
3301: - field - the field number
3303: Output Parameters:
3304: . sym - the symmetry describing the affect of orientation on the access of the data
3306: Level: developer
3308: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3309: @*/
3310: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3311: {
3314: 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);
3315: *sym = section->field[field]->sym;
3316: return(0);
3317: }
3319: /*@C
3320: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
3322: Not collective
3324: Input Parameters:
3325: + section - the section
3326: . numPoints - the number of points
3327: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3328: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3329: context, see DMPlexGetConeOrientation()).
3331: Output Parameter:
3332: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3333: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3334: identity).
3336: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3337: .vb
3338: const PetscInt **perms;
3339: const PetscScalar **rots;
3340: PetscInt lOffset;
3342: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3343: for (i = 0, lOffset = 0; i < numPoints; i++) {
3344: PetscInt point = points[2*i], dof, sOffset;
3345: const PetscInt *perm = perms ? perms[i] : NULL;
3346: const PetscScalar *rot = rots ? rots[i] : NULL;
3348: PetscSectionGetDof(section,point,&dof);
3349: PetscSectionGetOffset(section,point,&sOffset);
3351: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3352: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3353: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3354: lOffset += dof;
3355: }
3356: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3357: .ve
3359: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3360: .vb
3361: const PetscInt **perms;
3362: const PetscScalar **rots;
3363: PetscInt lOffset;
3365: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3366: for (i = 0, lOffset = 0; i < numPoints; i++) {
3367: PetscInt point = points[2*i], dof, sOffset;
3368: const PetscInt *perm = perms ? perms[i] : NULL;
3369: const PetscScalar *rot = rots ? rots[i] : NULL;
3371: PetscSectionGetDof(section,point,&dof);
3372: PetscSectionGetOffset(section,point,&sOff);
3374: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3375: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3376: offset += dof;
3377: }
3378: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3379: .ve
3381: Level: developer
3383: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3384: @*/
3385: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3386: {
3387: PetscSectionSym sym;
3388: PetscErrorCode ierr;
3393: if (perms) *perms = NULL;
3394: if (rots) *rots = NULL;
3395: sym = section->sym;
3396: if (sym && (perms || rots)) {
3397: SymWorkLink link;
3399: if (sym->workin) {
3400: link = sym->workin;
3401: sym->workin = sym->workin->next;
3402: } else {
3403: PetscNewLog(sym,&link);
3404: }
3405: if (numPoints > link->numPoints) {
3406: PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3407: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3408: link->numPoints = numPoints;
3409: }
3410: link->next = sym->workout;
3411: sym->workout = link;
3412: PetscArrayzero((PetscInt**)link->perms,numPoints);
3413: PetscArrayzero((PetscInt**)link->rots,numPoints);
3414: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3415: if (perms) *perms = link->perms;
3416: if (rots) *rots = link->rots;
3417: }
3418: return(0);
3419: }
3421: /*@C
3422: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
3424: Not collective
3426: Input Parameters:
3427: + section - the section
3428: . numPoints - the number of points
3429: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3430: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3431: context, see DMPlexGetConeOrientation()).
3433: Output Parameter:
3434: + perms - The permutations for the given orientations: set to NULL at conclusion
3435: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3437: Level: developer
3439: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3440: @*/
3441: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3442: {
3443: PetscSectionSym sym;
3447: sym = section->sym;
3448: if (sym && (perms || rots)) {
3449: SymWorkLink *p,link;
3451: for (p=&sym->workout; (link=*p); p=&link->next) {
3452: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3453: *p = link->next;
3454: link->next = sym->workin;
3455: sym->workin = link;
3456: if (perms) *perms = NULL;
3457: if (rots) *rots = NULL;
3458: return(0);
3459: }
3460: }
3461: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3462: }
3463: return(0);
3464: }
3466: /*@C
3467: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
3469: Not collective
3471: Input Parameters:
3472: + section - the section
3473: . field - the field of the section
3474: . numPoints - the number of points
3475: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3476: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3477: context, see DMPlexGetConeOrientation()).
3479: Output Parameter:
3480: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3481: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3482: identity).
3484: Level: developer
3486: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3487: @*/
3488: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3489: {
3494: 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);
3495: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3496: return(0);
3497: }
3499: /*@C
3500: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
3502: Not collective
3504: Input Parameters:
3505: + section - the section
3506: . field - the field number
3507: . numPoints - the number of points
3508: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3509: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3510: context, see DMPlexGetConeOrientation()).
3512: Output Parameter:
3513: + perms - The permutations for the given orientations: set to NULL at conclusion
3514: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3516: Level: developer
3518: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3519: @*/
3520: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3521: {
3526: 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);
3527: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3528: return(0);
3529: }
3531: /*@
3532: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3534: Not collective
3536: Input Parameter:
3537: . s - the global PetscSection
3539: Output Parameters:
3540: . flg - the flag
3542: Level: developer
3544: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3545: @*/
3546: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3547: {
3550: *flg = s->useFieldOff;
3551: return(0);
3552: }
3554: /*@
3555: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3557: Not collective
3559: Input Parameters:
3560: + s - the global PetscSection
3561: - flg - the flag
3563: Level: developer
3565: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3566: @*/
3567: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3568: {
3571: s->useFieldOff = flg;
3572: return(0);
3573: }
3575: #define PetscSectionExpandPoints_Loop(TYPE) \
3576: { \
3577: PetscInt i, n, o0, o1, size; \
3578: TYPE *a0 = (TYPE*)origArray, *a1; \
3579: PetscSectionGetStorageSize(s, &size); \
3580: PetscMalloc1(size, &a1); \
3581: for (i=0; i<npoints; i++) { \
3582: PetscSectionGetOffset(origSection, points_[i], &o0); \
3583: PetscSectionGetOffset(s, i, &o1); \
3584: PetscSectionGetDof(s, i, &n); \
3585: PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3586: } \
3587: *newArray = (void*)a1; \
3588: }
3590: /*@
3591: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3593: Not collective
3595: Input Parameters:
3596: + origSection - the PetscSection describing the layout of the array
3597: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3598: . origArray - the array; its size must be equal to the storage size of origSection
3599: - points - IS with points to extract; its indices must lie in the chart of origSection
3601: Output Parameters:
3602: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3603: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3605: Level: developer
3607: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3608: @*/
3609: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3610: {
3611: PetscSection s;
3612: const PetscInt *points_;
3613: PetscInt i, n, npoints, pStart, pEnd;
3614: PetscMPIInt unitsize;
3615: PetscErrorCode ierr;
3623: MPI_Type_size(dataType, &unitsize);
3624: ISGetLocalSize(points, &npoints);
3625: ISGetIndices(points, &points_);
3626: PetscSectionGetChart(origSection, &pStart, &pEnd);
3627: PetscSectionCreate(PETSC_COMM_SELF, &s);
3628: PetscSectionSetChart(s, 0, npoints);
3629: for (i=0; i<npoints; i++) {
3630: 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);
3631: PetscSectionGetDof(origSection, points_[i], &n);
3632: PetscSectionSetDof(s, i, n);
3633: }
3634: PetscSectionSetUp(s);
3635: if (newArray) {
3636: if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);}
3637: else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);}
3638: else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);}
3639: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3640: }
3641: if (newSection) {
3642: *newSection = s;
3643: } else {
3644: PetscSectionDestroy(&s);
3645: }
3646: return(0);
3647: }