Actual source code: section.c
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsf.h>
8: PetscClassId PETSC_SECTION_CLASSID;
10: /*@
11: PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
13: Collective
15: Input Parameters:
16: + comm - the MPI communicator
17: - s - pointer to the section
19: Level: beginner
21: Notes:
22: Typical calling sequence
23: $ PetscSectionCreate(MPI_Comm,PetscSection *);
24: $ PetscSectionSetNumFields(PetscSection, numFields);
25: $ PetscSectionSetChart(PetscSection,low,high);
26: $ PetscSectionSetDof(PetscSection,point,numdof);
27: $ PetscSectionSetUp(PetscSection);
28: $ PetscSectionGetOffset(PetscSection,point,PetscInt *);
29: $ PetscSectionDestroy(PetscSection);
31: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
32: recommended they not be used in user codes unless you really gain something in their use.
34: .seealso: PetscSection, PetscSectionDestroy()
35: @*/
36: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
37: {
42: ISInitializePackage();
44: PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
46: (*s)->pStart = -1;
47: (*s)->pEnd = -1;
48: (*s)->perm = NULL;
49: (*s)->pointMajor = PETSC_TRUE;
50: (*s)->includesConstraints = PETSC_TRUE;
51: (*s)->maxDof = 0;
52: (*s)->atlasDof = NULL;
53: (*s)->atlasOff = NULL;
54: (*s)->bc = NULL;
55: (*s)->bcIndices = NULL;
56: (*s)->setup = PETSC_FALSE;
57: (*s)->numFields = 0;
58: (*s)->fieldNames = NULL;
59: (*s)->field = NULL;
60: (*s)->useFieldOff = PETSC_FALSE;
61: (*s)->compNames = NULL;
62: (*s)->clObj = NULL;
63: (*s)->clHash = NULL;
64: (*s)->clSection = NULL;
65: (*s)->clPoints = NULL;
66: return(0);
67: }
69: /*@
70: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
72: Collective
74: Input Parameter:
75: . section - the PetscSection
77: Output Parameter:
78: . newSection - the copy
80: Level: intermediate
82: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
83: @*/
84: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
85: {
86: PetscSectionSym sym;
87: IS perm;
88: PetscInt numFields, f, c, pStart, pEnd, p;
89: PetscErrorCode ierr;
94: PetscSectionReset(newSection);
95: PetscSectionGetNumFields(section, &numFields);
96: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
97: for (f = 0; f < numFields; ++f) {
98: const char *fieldName = NULL, *compName = NULL;
99: PetscInt numComp = 0;
101: PetscSectionGetFieldName(section, f, &fieldName);
102: PetscSectionSetFieldName(newSection, f, fieldName);
103: PetscSectionGetFieldComponents(section, f, &numComp);
104: PetscSectionSetFieldComponents(newSection, f, numComp);
105: for (c = 0; c < numComp; ++c) {
106: PetscSectionGetComponentName(section, f, c, &compName);
107: PetscSectionSetComponentName(newSection, f, c, compName);
108: }
109: PetscSectionGetFieldSym(section, f, &sym);
110: PetscSectionSetFieldSym(newSection, f, sym);
111: }
112: PetscSectionGetChart(section, &pStart, &pEnd);
113: PetscSectionSetChart(newSection, pStart, pEnd);
114: PetscSectionGetPermutation(section, &perm);
115: PetscSectionSetPermutation(newSection, perm);
116: PetscSectionGetSym(section, &sym);
117: PetscSectionSetSym(newSection, sym);
118: for (p = pStart; p < pEnd; ++p) {
119: PetscInt dof, cdof, fcdof = 0;
121: PetscSectionGetDof(section, p, &dof);
122: PetscSectionSetDof(newSection, p, dof);
123: PetscSectionGetConstraintDof(section, p, &cdof);
124: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
125: for (f = 0; f < numFields; ++f) {
126: PetscSectionGetFieldDof(section, p, f, &dof);
127: PetscSectionSetFieldDof(newSection, p, f, dof);
128: if (cdof) {
129: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
130: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
131: }
132: }
133: }
134: PetscSectionSetUp(newSection);
135: for (p = pStart; p < pEnd; ++p) {
136: PetscInt off, cdof, fcdof = 0;
137: const PetscInt *cInd;
139: /* Must set offsets in case they do not agree with the prefix sums */
140: PetscSectionGetOffset(section, p, &off);
141: PetscSectionSetOffset(newSection, p, off);
142: PetscSectionGetConstraintDof(section, p, &cdof);
143: if (cdof) {
144: PetscSectionGetConstraintIndices(section, p, &cInd);
145: PetscSectionSetConstraintIndices(newSection, p, cInd);
146: for (f = 0; f < numFields; ++f) {
147: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
148: if (fcdof) {
149: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
150: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
151: }
152: }
153: }
154: }
155: return(0);
156: }
158: /*@
159: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
161: Collective
163: Input Parameter:
164: . section - the PetscSection
166: Output Parameter:
167: . newSection - the copy
169: Level: beginner
171: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
172: @*/
173: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
174: {
180: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
181: PetscSectionCopy(section, *newSection);
182: return(0);
183: }
185: /*@
186: PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database
188: Collective on PetscSection
190: Input Parameter:
191: . section - the PetscSection
193: Options Database:
194: . -petscsection_point_major the dof order
196: Level: intermediate
198: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
199: @*/
200: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
201: {
206: PetscObjectOptionsBegin((PetscObject) s);
207: PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
208: /* process any options handlers added with PetscObjectAddOptionsHandler() */
209: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
210: PetscOptionsEnd();
211: PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
212: return(0);
213: }
215: /*@
216: PetscSectionCompare - Compares two sections
218: Collective on PetscSection
220: Input Parameters:
221: + s1 - the first PetscSection
222: - s2 - the second PetscSection
224: Output Parameter:
225: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
227: Level: intermediate
229: Notes:
230: Field names are disregarded.
232: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
233: @*/
234: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
235: {
236: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
237: const PetscInt *idx1, *idx2;
238: IS perm1, perm2;
239: PetscBool flg;
240: PetscMPIInt mflg;
247: flg = PETSC_FALSE;
249: MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
250: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
251: *congruent = PETSC_FALSE;
252: return(0);
253: }
255: PetscSectionGetChart(s1, &pStart, &pEnd);
256: PetscSectionGetChart(s2, &n1, &n2);
257: if (pStart != n1 || pEnd != n2) goto not_congruent;
259: PetscSectionGetPermutation(s1, &perm1);
260: PetscSectionGetPermutation(s2, &perm2);
261: if (perm1 && perm2) {
262: ISEqual(perm1, perm2, congruent);
263: if (!(*congruent)) goto not_congruent;
264: } else if (perm1 != perm2) goto not_congruent;
266: for (p = pStart; p < pEnd; ++p) {
267: PetscSectionGetOffset(s1, p, &n1);
268: PetscSectionGetOffset(s2, p, &n2);
269: if (n1 != n2) goto not_congruent;
271: PetscSectionGetDof(s1, p, &n1);
272: PetscSectionGetDof(s2, p, &n2);
273: if (n1 != n2) goto not_congruent;
275: PetscSectionGetConstraintDof(s1, p, &ncdof);
276: PetscSectionGetConstraintDof(s2, p, &n2);
277: if (ncdof != n2) goto not_congruent;
279: PetscSectionGetConstraintIndices(s1, p, &idx1);
280: PetscSectionGetConstraintIndices(s2, p, &idx2);
281: PetscArraycmp(idx1, idx2, ncdof, congruent);
282: if (!(*congruent)) goto not_congruent;
283: }
285: PetscSectionGetNumFields(s1, &nfields);
286: PetscSectionGetNumFields(s2, &n2);
287: if (nfields != n2) goto not_congruent;
289: for (f = 0; f < nfields; ++f) {
290: PetscSectionGetFieldComponents(s1, f, &n1);
291: PetscSectionGetFieldComponents(s2, f, &n2);
292: if (n1 != n2) goto not_congruent;
294: for (p = pStart; p < pEnd; ++p) {
295: PetscSectionGetFieldOffset(s1, p, f, &n1);
296: PetscSectionGetFieldOffset(s2, p, f, &n2);
297: if (n1 != n2) goto not_congruent;
299: PetscSectionGetFieldDof(s1, p, f, &n1);
300: PetscSectionGetFieldDof(s2, p, f, &n2);
301: if (n1 != n2) goto not_congruent;
303: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
304: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
305: if (nfcdof != n2) goto not_congruent;
307: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
308: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
309: PetscArraycmp(idx1, idx2, nfcdof, congruent);
310: if (!(*congruent)) goto not_congruent;
311: }
312: }
314: flg = PETSC_TRUE;
315: not_congruent:
316: MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
317: return(0);
318: }
320: /*@
321: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
323: Not collective
325: Input Parameter:
326: . s - the PetscSection
328: Output Parameter:
329: . numFields - the number of fields defined, or 0 if none were defined
331: Level: intermediate
333: .seealso: PetscSectionSetNumFields()
334: @*/
335: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
336: {
340: *numFields = s->numFields;
341: return(0);
342: }
344: /*@
345: PetscSectionSetNumFields - Sets the number of fields.
347: Not collective
349: Input Parameters:
350: + s - the PetscSection
351: - numFields - the number of fields
353: Level: intermediate
355: .seealso: PetscSectionGetNumFields()
356: @*/
357: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
358: {
359: PetscInt f;
364: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %D must be positive", numFields);
365: PetscSectionReset(s);
367: s->numFields = numFields;
368: PetscMalloc1(s->numFields, &s->numFieldComponents);
369: PetscMalloc1(s->numFields, &s->fieldNames);
370: PetscMalloc1(s->numFields, &s->compNames);
371: PetscMalloc1(s->numFields, &s->field);
372: for (f = 0; f < s->numFields; ++f) {
373: char name[64];
375: s->numFieldComponents[f] = 1;
377: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
378: PetscSNPrintf(name, 64, "Field_%D", f);
379: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
380: PetscSNPrintf(name, 64, "Component_0");
381: PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
382: PetscStrallocpy(name, (char **) &s->compNames[f][0]);
383: }
384: return(0);
385: }
387: /*@C
388: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
390: Not Collective
392: Input Parameters:
393: + s - the PetscSection
394: - field - the field number
396: Output Parameter:
397: . fieldName - the field name
399: Level: intermediate
401: .seealso: PetscSectionSetFieldName()
402: @*/
403: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
404: {
408: 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);
409: *fieldName = s->fieldNames[field];
410: return(0);
411: }
413: /*@C
414: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
416: Not Collective
418: Input Parameters:
419: + s - the PetscSection
420: . field - the field number
421: - fieldName - the field name
423: Level: intermediate
425: .seealso: PetscSectionGetFieldName()
426: @*/
427: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
428: {
434: 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);
435: PetscFree(s->fieldNames[field]);
436: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
437: return(0);
438: }
440: /*@C
441: PetscSectionGetComponentName - Gets the name of a field component in the PetscSection
443: Not Collective
445: Input Parameters:
446: + s - the PetscSection
447: . field - the field number
448: . comp - the component number
449: - compName - the component name
451: Level: intermediate
453: .seealso: PetscSectionSetComponentName()
454: @*/
455: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
456: {
460: 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);
461: 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]);
462: *compName = s->compNames[field][comp];
463: return(0);
464: }
466: /*@C
467: PetscSectionSetComponentName - Sets the name of a field component in the PetscSection
469: Not Collective
471: Input Parameters:
472: + s - the PetscSection
473: . field - the field number
474: . comp - the component number
475: - compName - the component name
477: Level: intermediate
479: .seealso: PetscSectionGetComponentName()
480: @*/
481: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
482: {
488: 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);
489: 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]);
490: PetscFree(s->compNames[field][comp]);
491: PetscStrallocpy(compName, (char**) &s->compNames[field][comp]);
492: return(0);
493: }
495: /*@
496: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
498: Not collective
500: Input Parameters:
501: + s - the PetscSection
502: - field - the field number
504: Output Parameter:
505: . numComp - the number of field components
507: Level: intermediate
509: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
510: @*/
511: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
512: {
516: 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);
517: *numComp = s->numFieldComponents[field];
518: return(0);
519: }
521: /*@
522: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
524: Not collective
526: Input Parameters:
527: + s - the PetscSection
528: . field - the field number
529: - numComp - the number of field components
531: Level: intermediate
533: .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
534: @*/
535: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
536: {
538: PetscInt c;
539: char name[64];
543: 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);
544: if (s->compNames) {
545: for (c = 0; c < s->numFieldComponents[field]; ++c) {
546: PetscFree(s->compNames[field][c]);
547: }
548: PetscFree(s->compNames[field]);
549: }
551: s->numFieldComponents[field] = numComp;
552: if (numComp) {
553: PetscMalloc1(numComp, (char ***) &s->compNames[field]);
554: for (c = 0; c < numComp; ++c) {
555: PetscSNPrintf(name, 64, "%D", c);
556: PetscStrallocpy(name, (char **) &s->compNames[field][c]);
557: }
558: }
559: return(0);
560: }
562: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
563: {
567: if (!s->bc) {
568: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
569: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
570: }
571: return(0);
572: }
574: /*@
575: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
577: Not collective
579: Input Parameter:
580: . s - the PetscSection
582: Output Parameters:
583: + pStart - the first point
584: - pEnd - one past the last point
586: Level: intermediate
588: .seealso: PetscSectionSetChart(), PetscSectionCreate()
589: @*/
590: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
591: {
594: if (pStart) *pStart = s->pStart;
595: if (pEnd) *pEnd = s->pEnd;
596: return(0);
597: }
599: /*@
600: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
602: Not collective
604: Input Parameters:
605: + s - the PetscSection
606: . pStart - the first point
607: - pEnd - one past the last point
609: Level: intermediate
611: .seealso: PetscSectionGetChart(), PetscSectionCreate()
612: @*/
613: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
614: {
615: PetscInt f;
620: if (pStart == s->pStart && pEnd == s->pEnd) return(0);
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: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets
742: Not collective
744: Input Parameter:
745: . s - the PetscSection
747: Output Parameter:
748: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
750: Level: intermediate
752: .seealso: PetscSectionSetIncludesConstraints()
753: @*/
754: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
755: {
759: *includesConstraints = s->includesConstraints;
760: return(0);
761: }
763: /*@
764: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
766: Not collective
768: Input Parameters:
769: + s - the PetscSection
770: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
772: Not collective
774: Level: intermediate
776: .seealso: PetscSectionGetIncludesConstraints()
777: @*/
778: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
779: {
782: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
783: s->includesConstraints = includesConstraints;
784: return(0);
785: }
787: /*@
788: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
790: Not collective
792: Input Parameters:
793: + s - the PetscSection
794: - point - the point
796: Output Parameter:
797: . numDof - the number of dof
799: Level: intermediate
801: .seealso: PetscSectionSetDof(), PetscSectionCreate()
802: @*/
803: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
804: {
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: *numDof = s->atlasDof[point - s->pStart];
812: return(0);
813: }
815: /*@
816: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
818: Not collective
820: Input Parameters:
821: + s - the PetscSection
822: . point - the point
823: - numDof - the number of dof
825: Level: intermediate
827: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
828: @*/
829: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
830: {
833: 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);
834: s->atlasDof[point - s->pStart] = numDof;
835: return(0);
836: }
838: /*@
839: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
841: Not collective
843: Input Parameters:
844: + s - the PetscSection
845: . point - the point
846: - numDof - the number of additional dof
848: Level: intermediate
850: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
851: @*/
852: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
853: {
856: if (PetscDefined(USE_DEBUG)) {
857: 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);
858: }
859: s->atlasDof[point - s->pStart] += numDof;
860: return(0);
861: }
863: /*@
864: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
866: Not collective
868: Input Parameters:
869: + s - the PetscSection
870: . point - the point
871: - field - the field
873: Output Parameter:
874: . numDof - the number of dof
876: Level: intermediate
878: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
879: @*/
880: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
881: {
886: 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);
887: PetscSectionGetDof(s->field[field], point, numDof);
888: return(0);
889: }
891: /*@
892: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
894: Not collective
896: Input Parameters:
897: + s - the PetscSection
898: . point - the point
899: . field - the field
900: - numDof - the number of dof
902: Level: intermediate
904: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
905: @*/
906: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
907: {
912: 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);
913: PetscSectionSetDof(s->field[field], point, numDof);
914: return(0);
915: }
917: /*@
918: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
920: Not collective
922: Input Parameters:
923: + s - the PetscSection
924: . point - the point
925: . field - the field
926: - numDof - the number of dof
928: Level: intermediate
930: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
931: @*/
932: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
933: {
938: 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);
939: PetscSectionAddDof(s->field[field], point, numDof);
940: return(0);
941: }
943: /*@
944: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
946: Not collective
948: Input Parameters:
949: + s - the PetscSection
950: - point - the point
952: Output Parameter:
953: . numDof - the number of dof which are fixed by constraints
955: Level: intermediate
957: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
958: @*/
959: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
960: {
966: if (s->bc) {
967: PetscSectionGetDof(s->bc, point, numDof);
968: } else *numDof = 0;
969: return(0);
970: }
972: /*@
973: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
975: Not collective
977: Input Parameters:
978: + s - the PetscSection
979: . point - the point
980: - numDof - the number of dof which are fixed by constraints
982: Level: intermediate
984: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
985: @*/
986: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
987: {
992: if (numDof) {
993: PetscSectionCheckConstraints_Static(s);
994: PetscSectionSetDof(s->bc, point, numDof);
995: }
996: return(0);
997: }
999: /*@
1000: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
1002: Not collective
1004: Input Parameters:
1005: + s - the PetscSection
1006: . point - the point
1007: - numDof - the number of additional dof which are fixed by constraints
1009: Level: intermediate
1011: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
1012: @*/
1013: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1014: {
1019: if (numDof) {
1020: PetscSectionCheckConstraints_Static(s);
1021: PetscSectionAddDof(s->bc, point, numDof);
1022: }
1023: return(0);
1024: }
1026: /*@
1027: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
1029: Not collective
1031: Input Parameters:
1032: + s - the PetscSection
1033: . point - the point
1034: - field - the field
1036: Output Parameter:
1037: . numDof - the number of dof which are fixed by constraints
1039: Level: intermediate
1041: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
1042: @*/
1043: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1044: {
1050: 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);
1051: PetscSectionGetConstraintDof(s->field[field], point, numDof);
1052: return(0);
1053: }
1055: /*@
1056: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1058: Not collective
1060: Input Parameters:
1061: + s - the PetscSection
1062: . point - the point
1063: . field - the field
1064: - numDof - the number of dof which are fixed by constraints
1066: Level: intermediate
1068: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1069: @*/
1070: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1071: {
1076: 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);
1077: PetscSectionSetConstraintDof(s->field[field], point, numDof);
1078: return(0);
1079: }
1081: /*@
1082: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1084: Not collective
1086: Input Parameters:
1087: + s - the PetscSection
1088: . point - the point
1089: . field - the field
1090: - numDof - the number of additional dof which are fixed by constraints
1092: Level: intermediate
1094: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1095: @*/
1096: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1097: {
1102: 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);
1103: PetscSectionAddConstraintDof(s->field[field], point, numDof);
1104: return(0);
1105: }
1107: /*@
1108: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1110: Not collective
1112: Input Parameter:
1113: . s - the PetscSection
1115: Level: advanced
1117: .seealso: PetscSectionSetUp(), PetscSectionCreate()
1118: @*/
1119: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1120: {
1125: if (s->bc) {
1126: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
1128: PetscSectionSetUp(s->bc);
1129: PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1130: }
1131: return(0);
1132: }
1134: /*@
1135: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
1137: Not collective
1139: Input Parameter:
1140: . s - the PetscSection
1142: Level: intermediate
1144: .seealso: PetscSectionCreate()
1145: @*/
1146: PetscErrorCode PetscSectionSetUp(PetscSection s)
1147: {
1148: const PetscInt *pind = NULL;
1149: PetscInt offset = 0, foff, p, f;
1150: PetscErrorCode ierr;
1154: if (s->setup) return(0);
1155: s->setup = PETSC_TRUE;
1156: /* Set offsets and field offsets for all points */
1157: /* Assume that all fields have the same chart */
1158: if (!s->includesConstraints) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE");
1159: if (s->perm) {ISGetIndices(s->perm, &pind);}
1160: if (s->pointMajor) {
1161: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1162: const PetscInt q = pind ? pind[p] : p;
1164: /* Set point offset */
1165: s->atlasOff[q] = offset;
1166: offset += s->atlasDof[q];
1167: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
1168: /* Set field offset */
1169: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1170: PetscSection sf = s->field[f];
1172: sf->atlasOff[q] = foff;
1173: foff += sf->atlasDof[q];
1174: }
1175: }
1176: } else {
1177: /* Set field offsets for all points */
1178: for (f = 0; f < s->numFields; ++f) {
1179: PetscSection sf = s->field[f];
1181: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1182: const PetscInt q = pind ? pind[p] : p;
1184: sf->atlasOff[q] = offset;
1185: offset += sf->atlasDof[q];
1186: }
1187: }
1188: /* Disable point offsets since these are unused */
1189: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1190: s->atlasOff[p] = -1;
1191: s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1192: }
1193: }
1194: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1195: /* Setup BC sections */
1196: PetscSectionSetUpBC(s);
1197: for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1198: return(0);
1199: }
1201: /*@
1202: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1204: Not collective
1206: Input Parameters:
1207: . s - the PetscSection
1209: Output Parameter:
1210: . maxDof - the maximum dof
1212: Level: intermediate
1214: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1215: @*/
1216: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1217: {
1221: *maxDof = s->maxDof;
1222: return(0);
1223: }
1225: /*@
1226: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
1228: Not collective
1230: Input Parameter:
1231: . s - the PetscSection
1233: Output Parameter:
1234: . size - the size of an array which can hold all the dofs
1236: Level: intermediate
1238: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1239: @*/
1240: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1241: {
1242: PetscInt p, n = 0;
1247: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1248: *size = n;
1249: return(0);
1250: }
1252: /*@
1253: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1255: Not collective
1257: Input Parameter:
1258: . s - the PetscSection
1260: Output Parameter:
1261: . size - the size of an array which can hold all unconstrained dofs
1263: Level: intermediate
1265: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1266: @*/
1267: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1268: {
1269: PetscInt p, n = 0;
1274: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1275: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1276: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1277: }
1278: *size = n;
1279: return(0);
1280: }
1282: /*@
1283: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1284: the local section and an SF describing the section point overlap.
1286: Input Parameters:
1287: + s - The PetscSection for the local field layout
1288: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1289: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1290: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1292: Output Parameter:
1293: . gsection - The PetscSection for the global field layout
1295: Note: This gives negative sizes and offsets to points not owned by this process
1297: Level: intermediate
1299: .seealso: PetscSectionCreate()
1300: @*/
1301: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1302: {
1303: PetscSection gs;
1304: const PetscInt *pind = NULL;
1305: PetscInt *recv = NULL, *neg = NULL;
1306: PetscInt pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1307: PetscInt numFields, f, numComponents;
1308: PetscErrorCode ierr;
1316: if (!s->pointMajor) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "No support for field major ordering");
1317: PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1318: PetscSectionGetNumFields(s, &numFields);
1319: if (numFields > 0) {PetscSectionSetNumFields(gs, numFields);}
1320: PetscSectionGetChart(s, &pStart, &pEnd);
1321: PetscSectionSetChart(gs, pStart, pEnd);
1322: gs->includesConstraints = includeConstraints;
1323: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1324: nlocal = nroots; /* The local/leaf space matches global/root space */
1325: /* Must allocate for all points visible to SF, which may be more than this section */
1326: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1327: PetscSFGetLeafRange(sf, NULL, &maxleaf);
1328: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1329: if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1330: PetscMalloc2(nroots,&neg,nlocal,&recv);
1331: PetscArrayzero(neg,nroots);
1332: }
1333: /* Mark all local points with negative dof */
1334: for (p = pStart; p < pEnd; ++p) {
1335: PetscSectionGetDof(s, p, &dof);
1336: PetscSectionSetDof(gs, p, dof);
1337: PetscSectionGetConstraintDof(s, p, &cdof);
1338: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1339: if (neg) neg[p] = -(dof+1);
1340: }
1341: PetscSectionSetUpBC(gs);
1342: 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]);}
1343: if (nroots >= 0) {
1344: PetscArrayzero(recv,nlocal);
1345: PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1346: PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1347: for (p = pStart; p < pEnd; ++p) {
1348: if (recv[p] < 0) {
1349: gs->atlasDof[p-pStart] = recv[p];
1350: PetscSectionGetDof(s, p, &dof);
1351: 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);
1352: }
1353: }
1354: }
1355: /* Calculate new sizes, get process offset, and calculate point offsets */
1356: if (s->perm) {ISGetIndices(s->perm, &pind);}
1357: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1358: const PetscInt q = pind ? pind[p] : p;
1360: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1361: gs->atlasOff[q] = off;
1362: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1363: }
1364: if (!localOffsets) {
1365: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1366: globalOff -= off;
1367: }
1368: for (p = pStart, off = 0; p < pEnd; ++p) {
1369: gs->atlasOff[p-pStart] += globalOff;
1370: if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1371: }
1372: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1373: /* Put in negative offsets for ghost points */
1374: if (nroots >= 0) {
1375: PetscArrayzero(recv,nlocal);
1376: PetscSFBcastBegin(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1377: PetscSFBcastEnd(sf, MPIU_INT, neg, recv,MPI_REPLACE);
1378: for (p = pStart; p < pEnd; ++p) {
1379: if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1380: }
1381: }
1382: PetscFree2(neg,recv);
1383: /* Set field dofs/offsets/constraints */
1384: for (f = 0; f < numFields; ++f) {
1385: gs->field[f]->includesConstraints = includeConstraints;
1386: PetscSectionGetFieldComponents(s, f, &numComponents);
1387: PetscSectionSetFieldComponents(gs, f, numComponents);
1388: }
1389: for (p = pStart; p < pEnd; ++p) {
1390: PetscSectionGetOffset(gs, p, &off);
1391: for (f = 0, foff = off; f < numFields; ++f) {
1392: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
1393: if (!includeConstraints && cdof > 0) {PetscSectionSetFieldConstraintDof(gs, p, f, cdof);}
1394: PetscSectionGetFieldDof(s, p, f, &dof);
1395: PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof);
1396: PetscSectionSetFieldOffset(gs, p, f, foff);
1397: PetscSectionGetFieldConstraintDof(gs, p, f, &cdof);
1398: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1399: }
1400: }
1401: for (f = 0; f < numFields; ++f) {
1402: PetscSection gfs = gs->field[f];
1404: PetscSectionSetUpBC(gfs);
1405: if (gfs->bcIndices) {PetscArraycpy(gfs->bcIndices, s->field[f]->bcIndices, gfs->bc->atlasOff[gfs->bc->pEnd-gfs->bc->pStart-1] + gfs->bc->atlasDof[gfs->bc->pEnd-gfs->bc->pStart-1]);}
1406: }
1407: gs->setup = PETSC_TRUE;
1408: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1409: *gsection = gs;
1410: return(0);
1411: }
1413: /*@
1414: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1415: the local section and an SF describing the section point overlap.
1417: Input Parameters:
1418: + s - The PetscSection for the local field layout
1419: . sf - The SF describing parallel layout of the section points
1420: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1421: . numExcludes - The number of exclusion ranges
1422: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1424: Output Parameter:
1425: . gsection - The PetscSection for the global field layout
1427: Note: This gives negative sizes and offsets to points not owned by this process
1429: Level: advanced
1431: .seealso: PetscSectionCreate()
1432: @*/
1433: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1434: {
1435: const PetscInt *pind = NULL;
1436: PetscInt *neg = NULL, *tmpOff = NULL;
1437: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1438: PetscErrorCode ierr;
1444: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1445: PetscSectionGetChart(s, &pStart, &pEnd);
1446: PetscSectionSetChart(*gsection, pStart, pEnd);
1447: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1448: if (nroots >= 0) {
1449: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1450: PetscCalloc1(nroots, &neg);
1451: if (nroots > pEnd-pStart) {
1452: PetscCalloc1(nroots, &tmpOff);
1453: } else {
1454: tmpOff = &(*gsection)->atlasDof[-pStart];
1455: }
1456: }
1457: /* Mark ghost points with negative dof */
1458: for (p = pStart; p < pEnd; ++p) {
1459: for (e = 0; e < numExcludes; ++e) {
1460: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1461: PetscSectionSetDof(*gsection, p, 0);
1462: break;
1463: }
1464: }
1465: if (e < numExcludes) continue;
1466: PetscSectionGetDof(s, p, &dof);
1467: PetscSectionSetDof(*gsection, p, dof);
1468: PetscSectionGetConstraintDof(s, p, &cdof);
1469: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1470: if (neg) neg[p] = -(dof+1);
1471: }
1472: PetscSectionSetUpBC(*gsection);
1473: if (nroots >= 0) {
1474: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1475: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1476: if (nroots > pEnd - pStart) {
1477: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1478: }
1479: }
1480: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1481: if (s->perm) {ISGetIndices(s->perm, &pind);}
1482: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1483: const PetscInt q = pind ? pind[p] : p;
1485: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1486: (*gsection)->atlasOff[q] = off;
1487: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1488: }
1489: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1490: globalOff -= off;
1491: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1492: (*gsection)->atlasOff[p] += globalOff;
1493: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1494: }
1495: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1496: /* Put in negative offsets for ghost points */
1497: if (nroots >= 0) {
1498: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1499: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1500: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff,MPI_REPLACE);
1501: if (nroots > pEnd - pStart) {
1502: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1503: }
1504: }
1505: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1506: PetscFree(neg);
1507: return(0);
1508: }
1510: /*@
1511: PetscSectionGetPointLayout - Get the PetscLayout associated with the section points
1513: Collective on comm
1515: Input Parameters:
1516: + comm - The MPI_Comm
1517: - s - The PetscSection
1519: Output Parameter:
1520: . layout - The point layout for the section
1522: Note: This is usually called for the default global section.
1524: Level: advanced
1526: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1527: @*/
1528: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1529: {
1530: PetscInt pStart, pEnd, p, localSize = 0;
1534: PetscSectionGetChart(s, &pStart, &pEnd);
1535: for (p = pStart; p < pEnd; ++p) {
1536: PetscInt dof;
1538: PetscSectionGetDof(s, p, &dof);
1539: if (dof > 0) ++localSize;
1540: }
1541: PetscLayoutCreate(comm, layout);
1542: PetscLayoutSetLocalSize(*layout, localSize);
1543: PetscLayoutSetBlockSize(*layout, 1);
1544: PetscLayoutSetUp(*layout);
1545: return(0);
1546: }
1548: /*@
1549: PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.
1551: Collective on comm
1553: Input Parameters:
1554: + comm - The MPI_Comm
1555: - s - The PetscSection
1557: Output Parameter:
1558: . layout - The dof layout for the section
1560: Note: This is usually called for the default global section.
1562: Level: advanced
1564: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1565: @*/
1566: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1567: {
1568: PetscInt pStart, pEnd, p, localSize = 0;
1574: PetscSectionGetChart(s, &pStart, &pEnd);
1575: for (p = pStart; p < pEnd; ++p) {
1576: PetscInt dof,cdof;
1578: PetscSectionGetDof(s, p, &dof);
1579: PetscSectionGetConstraintDof(s, p, &cdof);
1580: if (dof-cdof > 0) localSize += dof-cdof;
1581: }
1582: PetscLayoutCreate(comm, layout);
1583: PetscLayoutSetLocalSize(*layout, localSize);
1584: PetscLayoutSetBlockSize(*layout, 1);
1585: PetscLayoutSetUp(*layout);
1586: return(0);
1587: }
1589: /*@
1590: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1592: Not collective
1594: Input Parameters:
1595: + s - the PetscSection
1596: - point - the point
1598: Output Parameter:
1599: . offset - the offset
1601: Level: intermediate
1603: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1604: @*/
1605: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1606: {
1610: if (PetscDefined(USE_DEBUG)) {
1611: 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);
1612: }
1613: *offset = s->atlasOff[point - s->pStart];
1614: return(0);
1615: }
1617: /*@
1618: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1620: Not collective
1622: Input Parameters:
1623: + s - the PetscSection
1624: . point - the point
1625: - offset - the offset
1627: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1629: Level: intermediate
1631: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1632: @*/
1633: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1634: {
1637: 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);
1638: s->atlasOff[point - s->pStart] = offset;
1639: return(0);
1640: }
1642: /*@
1643: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1645: Not collective
1647: Input Parameters:
1648: + s - the PetscSection
1649: . point - the point
1650: - field - the field
1652: Output Parameter:
1653: . offset - the offset
1655: Level: intermediate
1657: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1658: @*/
1659: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1660: {
1666: 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);
1667: PetscSectionGetOffset(s->field[field], point, offset);
1668: return(0);
1669: }
1671: /*@
1672: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given field at a point.
1674: Not collective
1676: Input Parameters:
1677: + s - the PetscSection
1678: . point - the point
1679: . field - the field
1680: - offset - the offset
1682: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1684: Level: intermediate
1686: .seealso: PetscSectionGetFieldOffset(), PetscSectionSetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1687: @*/
1688: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1689: {
1694: 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);
1695: PetscSectionSetOffset(s->field[field], point, offset);
1696: return(0);
1697: }
1699: /*@
1700: PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.
1702: Not collective
1704: Input Parameters:
1705: + s - the PetscSection
1706: . point - the point
1707: - field - the field
1709: Output Parameter:
1710: . offset - the offset
1712: Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1713: this point, what number is the first dof with this field.
1715: Level: advanced
1717: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1718: @*/
1719: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1720: {
1721: PetscInt off, foff;
1727: 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);
1728: PetscSectionGetOffset(s, point, &off);
1729: PetscSectionGetOffset(s->field[field], point, &foff);
1730: *offset = foff - off;
1731: return(0);
1732: }
1734: /*@
1735: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1737: Not collective
1739: Input Parameter:
1740: . s - the PetscSection
1742: Output Parameters:
1743: + start - the minimum offset
1744: - end - one more than the maximum offset
1746: Level: intermediate
1748: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1749: @*/
1750: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1751: {
1752: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1757: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1758: PetscSectionGetChart(s, &pStart, &pEnd);
1759: for (p = 0; p < pEnd-pStart; ++p) {
1760: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1762: if (off >= 0) {
1763: os = PetscMin(os, off);
1764: oe = PetscMax(oe, off+dof);
1765: }
1766: }
1767: if (start) *start = os;
1768: if (end) *end = oe;
1769: return(0);
1770: }
1772: /*@
1773: PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields
1775: Collective on s
1777: Input Parameters:
1778: + s - the PetscSection
1779: . len - the number of subfields
1780: - fields - the subfield numbers
1782: Output Parameter:
1783: . subs - the subsection
1785: Note: The section offsets now refer to a new, smaller vector.
1787: Level: advanced
1789: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1790: @*/
1791: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1792: {
1793: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
1797: if (!len) return(0);
1801: PetscSectionGetNumFields(s, &nF);
1802: if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1803: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1804: PetscSectionSetNumFields(*subs, len);
1805: for (f = 0; f < len; ++f) {
1806: const char *name = NULL;
1807: PetscInt numComp = 0;
1809: PetscSectionGetFieldName(s, fields[f], &name);
1810: PetscSectionSetFieldName(*subs, f, name);
1811: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1812: PetscSectionSetFieldComponents(*subs, f, numComp);
1813: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1814: PetscSectionGetComponentName(s, fields[f], c, &name);
1815: PetscSectionSetComponentName(*subs, f, c, name);
1816: }
1817: }
1818: PetscSectionGetChart(s, &pStart, &pEnd);
1819: PetscSectionSetChart(*subs, pStart, pEnd);
1820: for (p = pStart; p < pEnd; ++p) {
1821: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1823: for (f = 0; f < len; ++f) {
1824: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1825: PetscSectionSetFieldDof(*subs, p, f, fdof);
1826: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1827: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1828: dof += fdof;
1829: cdof += cfdof;
1830: }
1831: PetscSectionSetDof(*subs, p, dof);
1832: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1833: maxCdof = PetscMax(cdof, maxCdof);
1834: }
1835: PetscSectionSetUp(*subs);
1836: if (maxCdof) {
1837: PetscInt *indices;
1839: PetscMalloc1(maxCdof, &indices);
1840: for (p = pStart; p < pEnd; ++p) {
1841: PetscInt cdof;
1843: PetscSectionGetConstraintDof(*subs, p, &cdof);
1844: if (cdof) {
1845: const PetscInt *oldIndices = NULL;
1846: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1848: for (f = 0; f < len; ++f) {
1849: PetscInt oldFoff = 0, oldf;
1851: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1852: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1853: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1854: /* This can be sped up if we assume sorted fields */
1855: for (oldf = 0; oldf < fields[f]; ++oldf) {
1856: PetscInt oldfdof = 0;
1857: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1858: oldFoff += oldfdof;
1859: }
1860: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1861: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1862: numConst += cfdof;
1863: fOff += fdof;
1864: }
1865: PetscSectionSetConstraintIndices(*subs, p, indices);
1866: }
1867: }
1868: PetscFree(indices);
1869: }
1870: return(0);
1871: }
1873: /*@
1874: PetscSectionCreateSupersection - Create a new, larger section composed of the input sections
1876: Collective on s
1878: Input Parameters:
1879: + s - the input sections
1880: - len - the number of input sections
1882: Output Parameter:
1883: . supers - the supersection
1885: Note: The section offsets now refer to a new, larger vector.
1887: Level: advanced
1889: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1890: @*/
1891: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1892: {
1893: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1897: if (!len) return(0);
1898: for (i = 0; i < len; ++i) {
1899: PetscInt nf, pStarti, pEndi;
1901: PetscSectionGetNumFields(s[i], &nf);
1902: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1903: pStart = PetscMin(pStart, pStarti);
1904: pEnd = PetscMax(pEnd, pEndi);
1905: Nf += nf;
1906: }
1907: PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1908: PetscSectionSetNumFields(*supers, Nf);
1909: for (i = 0, f = 0; i < len; ++i) {
1910: PetscInt nf, fi, ci;
1912: PetscSectionGetNumFields(s[i], &nf);
1913: for (fi = 0; fi < nf; ++fi, ++f) {
1914: const char *name = NULL;
1915: PetscInt numComp = 0;
1917: PetscSectionGetFieldName(s[i], fi, &name);
1918: PetscSectionSetFieldName(*supers, f, name);
1919: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1920: PetscSectionSetFieldComponents(*supers, f, numComp);
1921: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1922: PetscSectionGetComponentName(s[i], fi, ci, &name);
1923: PetscSectionSetComponentName(*supers, f, ci, name);
1924: }
1925: }
1926: }
1927: PetscSectionSetChart(*supers, pStart, pEnd);
1928: for (p = pStart; p < pEnd; ++p) {
1929: PetscInt dof = 0, cdof = 0;
1931: for (i = 0, f = 0; i < len; ++i) {
1932: PetscInt nf, fi, pStarti, pEndi;
1933: PetscInt fdof = 0, cfdof = 0;
1935: PetscSectionGetNumFields(s[i], &nf);
1936: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1937: if ((p < pStarti) || (p >= pEndi)) continue;
1938: for (fi = 0; fi < nf; ++fi, ++f) {
1939: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1940: PetscSectionAddFieldDof(*supers, p, f, fdof);
1941: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1942: if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1943: dof += fdof;
1944: cdof += cfdof;
1945: }
1946: }
1947: PetscSectionSetDof(*supers, p, dof);
1948: if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1949: maxCdof = PetscMax(cdof, maxCdof);
1950: }
1951: PetscSectionSetUp(*supers);
1952: if (maxCdof) {
1953: PetscInt *indices;
1955: PetscMalloc1(maxCdof, &indices);
1956: for (p = pStart; p < pEnd; ++p) {
1957: PetscInt cdof;
1959: PetscSectionGetConstraintDof(*supers, p, &cdof);
1960: if (cdof) {
1961: PetscInt dof, numConst = 0, fOff = 0;
1963: for (i = 0, f = 0; i < len; ++i) {
1964: const PetscInt *oldIndices = NULL;
1965: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1967: PetscSectionGetNumFields(s[i], &nf);
1968: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1969: if ((p < pStarti) || (p >= pEndi)) continue;
1970: for (fi = 0; fi < nf; ++fi, ++f) {
1971: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1972: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1973: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1974: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc];
1975: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1976: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] += fOff;
1977: numConst += cfdof;
1978: }
1979: PetscSectionGetDof(s[i], p, &dof);
1980: fOff += dof;
1981: }
1982: PetscSectionSetConstraintIndices(*supers, p, indices);
1983: }
1984: }
1985: PetscFree(indices);
1986: }
1987: return(0);
1988: }
1990: /*@
1991: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
1993: Collective on s
1995: Input Parameters:
1996: + s - the PetscSection
1997: - subpointMap - a sorted list of points in the original mesh which are in the submesh
1999: Output Parameter:
2000: . subs - the subsection
2002: Note: The section offsets now refer to a new, smaller vector.
2004: Level: advanced
2006: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
2007: @*/
2008: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2009: {
2010: const PetscInt *points = NULL, *indices = NULL;
2011: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;
2018: PetscSectionGetNumFields(s, &numFields);
2019: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
2020: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
2021: for (f = 0; f < numFields; ++f) {
2022: const char *name = NULL;
2023: PetscInt numComp = 0;
2025: PetscSectionGetFieldName(s, f, &name);
2026: PetscSectionSetFieldName(*subs, f, name);
2027: PetscSectionGetFieldComponents(s, f, &numComp);
2028: PetscSectionSetFieldComponents(*subs, f, numComp);
2029: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2030: PetscSectionGetComponentName(s, f, c, &name);
2031: PetscSectionSetComponentName(*subs, f, c, name);
2032: }
2033: }
2034: /* For right now, we do not try to squeeze the subchart */
2035: if (subpointMap) {
2036: ISGetSize(subpointMap, &numSubpoints);
2037: ISGetIndices(subpointMap, &points);
2038: }
2039: PetscSectionGetChart(s, &pStart, &pEnd);
2040: PetscSectionSetChart(*subs, 0, numSubpoints);
2041: for (p = pStart; p < pEnd; ++p) {
2042: PetscInt dof, cdof, fdof = 0, cfdof = 0;
2044: PetscFindInt(p, numSubpoints, points, &subp);
2045: if (subp < 0) continue;
2046: for (f = 0; f < numFields; ++f) {
2047: PetscSectionGetFieldDof(s, p, f, &fdof);
2048: PetscSectionSetFieldDof(*subs, subp, f, fdof);
2049: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
2050: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
2051: }
2052: PetscSectionGetDof(s, p, &dof);
2053: PetscSectionSetDof(*subs, subp, dof);
2054: PetscSectionGetConstraintDof(s, p, &cdof);
2055: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
2056: }
2057: PetscSectionSetUp(*subs);
2058: /* Change offsets to original offsets */
2059: for (p = pStart; p < pEnd; ++p) {
2060: PetscInt off, foff = 0;
2062: PetscFindInt(p, numSubpoints, points, &subp);
2063: if (subp < 0) continue;
2064: for (f = 0; f < numFields; ++f) {
2065: PetscSectionGetFieldOffset(s, p, f, &foff);
2066: PetscSectionSetFieldOffset(*subs, subp, f, foff);
2067: }
2068: PetscSectionGetOffset(s, p, &off);
2069: PetscSectionSetOffset(*subs, subp, off);
2070: }
2071: /* Copy constraint indices */
2072: for (subp = 0; subp < numSubpoints; ++subp) {
2073: PetscInt cdof;
2075: PetscSectionGetConstraintDof(*subs, subp, &cdof);
2076: if (cdof) {
2077: for (f = 0; f < numFields; ++f) {
2078: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
2079: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2080: }
2081: PetscSectionGetConstraintIndices(s, points[subp], &indices);
2082: PetscSectionSetConstraintIndices(*subs, subp, indices);
2083: }
2084: }
2085: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
2086: return(0);
2087: }
2089: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2090: {
2091: PetscInt p;
2092: PetscMPIInt rank;
2096: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2097: PetscViewerASCIIPushSynchronized(viewer);
2098: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2099: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2100: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2101: PetscInt b;
2103: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2104: if (s->bcIndices) {
2105: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
2106: PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
2107: }
2108: }
2109: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2110: } else {
2111: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2112: }
2113: }
2114: PetscViewerFlush(viewer);
2115: PetscViewerASCIIPopSynchronized(viewer);
2116: if (s->sym) {
2117: PetscViewerASCIIPushTab(viewer);
2118: PetscSectionSymView(s->sym,viewer);
2119: PetscViewerASCIIPopTab(viewer);
2120: }
2121: return(0);
2122: }
2124: /*@C
2125: PetscSectionViewFromOptions - View from Options
2127: Collective on PetscSection
2129: Input Parameters:
2130: + A - the PetscSection object to view
2131: . obj - Optional object
2132: - name - command line option
2134: Level: intermediate
2135: .seealso: PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
2136: @*/
2137: PetscErrorCode PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
2138: {
2143: PetscObjectViewFromOptions((PetscObject)A,obj,name);
2144: return(0);
2145: }
2147: /*@C
2148: PetscSectionView - Views a PetscSection
2150: Collective on PetscSection
2152: Input Parameters:
2153: + s - the PetscSection object to view
2154: - v - the viewer
2156: Note:
2157: PetscSectionView(), when viewer is of type PETSCVIEWERHDF5, only saves
2158: distribution independent data, such as dofs, offsets, constraint dofs,
2159: and constraint indices. Points that have negative dofs, for instance,
2160: are not saved as they represent points owned by other processes.
2161: Point numbering and rank assignment is currently not stored.
2162: The saved section can be loaded with PetscSectionLoad().
2164: Level: beginner
2166: .seealso PetscSectionCreate(), PetscSectionDestroy(), PetscSectionLoad()
2167: @*/
2168: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2169: {
2170: PetscBool isascii, ishdf5;
2171: PetscInt f;
2176: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
2178: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
2179: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
2180: if (isascii) {
2181: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
2182: if (s->numFields) {
2183: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
2184: for (f = 0; f < s->numFields; ++f) {
2185: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
2186: PetscSectionView_ASCII(s->field[f], viewer);
2187: }
2188: } else {
2189: PetscSectionView_ASCII(s, viewer);
2190: }
2191: } else if (ishdf5) {
2192: #if PetscDefined(HAVE_HDF5)
2193: PetscSectionView_HDF5_Internal(s, viewer);
2194: #else
2195: SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2196: #endif
2197: }
2198: return(0);
2199: }
2201: /*@C
2202: PetscSectionLoad - Loads a PetscSection
2204: Collective on PetscSection
2206: Input Parameters:
2207: + s - the PetscSection object to load
2208: - v - the viewer
2210: Note:
2211: PetscSectionLoad(), when viewer is of type PETSCVIEWERHDF5, loads
2212: a section saved with PetscSectionView(). The number of processes
2213: used here (N) does not need to be the same as that used when saving.
2214: After calling this function, the chart of s on rank i will be set
2215: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2216: saved section points.
2218: Level: beginner
2220: .seealso PetscSectionCreate(), PetscSectionDestroy(), PetscSectionView()
2221: @*/
2222: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2223: {
2224: PetscBool ishdf5;
2230: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
2231: if (ishdf5) {
2232: #if PetscDefined(HAVE_HDF5)
2233: PetscSectionLoad_HDF5_Internal(s, viewer);
2234: return(0);
2235: #else
2236: SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2237: #endif
2238: } else SETERRQ1(PetscObjectComm((PetscObject) s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2239: }
2241: static PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2242: {
2244: PetscSectionClosurePermVal clVal;
2247: if (!section->clHash) return(0);
2248: kh_foreach_value(section->clHash, clVal, {
2249: PetscFree(clVal.perm);
2250: PetscFree(clVal.invPerm);
2251: });
2252: kh_destroy(ClPerm, section->clHash);
2253: section->clHash = NULL;
2254: return(0);
2255: }
2257: /*@
2258: PetscSectionReset - Frees all section data.
2260: Not collective
2262: Input Parameters:
2263: . s - the PetscSection
2265: Level: beginner
2267: .seealso: PetscSection, PetscSectionCreate()
2268: @*/
2269: PetscErrorCode PetscSectionReset(PetscSection s)
2270: {
2271: PetscInt f, c;
2276: for (f = 0; f < s->numFields; ++f) {
2277: PetscSectionDestroy(&s->field[f]);
2278: PetscFree(s->fieldNames[f]);
2279: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2280: PetscFree(s->compNames[f][c]);
2281: }
2282: PetscFree(s->compNames[f]);
2283: }
2284: PetscFree(s->numFieldComponents);
2285: PetscFree(s->fieldNames);
2286: PetscFree(s->compNames);
2287: PetscFree(s->field);
2288: PetscSectionDestroy(&s->bc);
2289: PetscFree(s->bcIndices);
2290: PetscFree2(s->atlasDof, s->atlasOff);
2291: PetscSectionDestroy(&s->clSection);
2292: ISDestroy(&s->clPoints);
2293: ISDestroy(&s->perm);
2294: PetscSectionResetClosurePermutation(s);
2295: PetscSectionSymDestroy(&s->sym);
2296: PetscSectionDestroy(&s->clSection);
2297: ISDestroy(&s->clPoints);
2299: s->pStart = -1;
2300: s->pEnd = -1;
2301: s->maxDof = 0;
2302: s->setup = PETSC_FALSE;
2303: s->numFields = 0;
2304: s->clObj = NULL;
2305: return(0);
2306: }
2308: /*@
2309: PetscSectionDestroy - Frees a section object and frees its range if that exists.
2311: Not collective
2313: Input Parameters:
2314: . s - the PetscSection
2316: Level: beginner
2318: .seealso: PetscSection, PetscSectionCreate()
2319: @*/
2320: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2321: {
2325: if (!*s) return(0);
2327: if (--((PetscObject)(*s))->refct > 0) {
2328: *s = NULL;
2329: return(0);
2330: }
2331: PetscSectionReset(*s);
2332: PetscHeaderDestroy(s);
2333: return(0);
2334: }
2336: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2337: {
2338: const PetscInt p = point - s->pStart;
2342: *values = &baseArray[s->atlasOff[p]];
2343: return(0);
2344: }
2346: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2347: {
2348: PetscInt *array;
2349: const PetscInt p = point - s->pStart;
2350: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2351: PetscInt cDim = 0;
2356: PetscSectionGetConstraintDof(s, p, &cDim);
2357: array = &baseArray[s->atlasOff[p]];
2358: if (!cDim) {
2359: if (orientation >= 0) {
2360: const PetscInt dim = s->atlasDof[p];
2361: PetscInt i;
2363: if (mode == INSERT_VALUES) {
2364: for (i = 0; i < dim; ++i) array[i] = values[i];
2365: } else {
2366: for (i = 0; i < dim; ++i) array[i] += values[i];
2367: }
2368: } else {
2369: PetscInt offset = 0;
2370: PetscInt j = -1, field, i;
2372: for (field = 0; field < s->numFields; ++field) {
2373: const PetscInt dim = s->field[field]->atlasDof[p];
2375: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2376: offset += dim;
2377: }
2378: }
2379: } else {
2380: if (orientation >= 0) {
2381: const PetscInt dim = s->atlasDof[p];
2382: PetscInt cInd = 0, i;
2383: const PetscInt *cDof;
2385: PetscSectionGetConstraintIndices(s, point, &cDof);
2386: if (mode == INSERT_VALUES) {
2387: for (i = 0; i < dim; ++i) {
2388: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2389: array[i] = values[i];
2390: }
2391: } else {
2392: for (i = 0; i < dim; ++i) {
2393: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2394: array[i] += values[i];
2395: }
2396: }
2397: } else {
2398: const PetscInt *cDof;
2399: PetscInt offset = 0;
2400: PetscInt cOffset = 0;
2401: PetscInt j = 0, field;
2403: PetscSectionGetConstraintIndices(s, point, &cDof);
2404: for (field = 0; field < s->numFields; ++field) {
2405: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2406: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2407: const PetscInt sDim = dim - tDim;
2408: PetscInt cInd = 0, i,k;
2410: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2411: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2412: array[j] = values[k];
2413: }
2414: offset += dim;
2415: cOffset += dim - tDim;
2416: }
2417: }
2418: }
2419: return(0);
2420: }
2422: /*@C
2423: PetscSectionHasConstraints - Determine whether a section has constrained dofs
2425: Not collective
2427: Input Parameter:
2428: . s - The PetscSection
2430: Output Parameter:
2431: . hasConstraints - flag indicating that the section has constrained dofs
2433: Level: intermediate
2435: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2436: @*/
2437: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2438: {
2442: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2443: return(0);
2444: }
2446: /*@C
2447: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2449: Not collective
2451: Input Parameters:
2452: + s - The PetscSection
2453: - point - The point
2455: Output Parameter:
2456: . indices - The constrained dofs
2458: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
2460: Level: intermediate
2462: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2463: @*/
2464: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2465: {
2470: if (s->bc) {
2471: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2472: } else *indices = NULL;
2473: return(0);
2474: }
2476: /*@C
2477: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2479: Not collective
2481: Input Parameters:
2482: + s - The PetscSection
2483: . point - The point
2484: - indices - The constrained dofs
2486: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
2488: Level: intermediate
2490: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2491: @*/
2492: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2493: {
2498: if (s->bc) {
2499: const PetscInt dof = s->atlasDof[point];
2500: const PetscInt cdof = s->bc->atlasDof[point];
2501: PetscInt d;
2503: for (d = 0; d < cdof; ++d) {
2504: if (indices[d] >= dof) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D dof %D, invalid constraint index[%D]: %D", point, dof, d, indices[d]);
2505: }
2506: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2507: }
2508: return(0);
2509: }
2511: /*@C
2512: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2514: Not collective
2516: Input Parameters:
2517: + s - The PetscSection
2518: . field - The field number
2519: - point - The point
2521: Output Parameter:
2522: . indices - The constrained dofs sorted in ascending order
2524: Notes:
2525: 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().
2527: Fortran Note:
2528: In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()
2530: Level: intermediate
2532: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2533: @*/
2534: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2535: {
2540: 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);
2541: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2542: return(0);
2543: }
2545: /*@C
2546: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2548: Not collective
2550: Input Parameters:
2551: + s - The PetscSection
2552: . point - The point
2553: . field - The field number
2554: - indices - The constrained dofs
2556: Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()
2558: Level: intermediate
2560: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2561: @*/
2562: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2563: {
2568: 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);
2569: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2570: return(0);
2571: }
2573: /*@
2574: PetscSectionPermute - Reorder the section according to the input point permutation
2576: Collective on PetscSection
2578: Input Parameters:
2579: + section - The PetscSection object
2580: - perm - The point permutation, old point p becomes new point perm[p]
2582: Output Parameter:
2583: . sectionNew - The permuted PetscSection
2585: Level: intermediate
2587: .seealso: MatPermute()
2588: @*/
2589: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2590: {
2591: PetscSection s = section, sNew;
2592: const PetscInt *perm;
2593: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2594: PetscErrorCode ierr;
2600: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2601: PetscSectionGetNumFields(s, &numFields);
2602: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2603: for (f = 0; f < numFields; ++f) {
2604: const char *name;
2605: PetscInt numComp;
2607: PetscSectionGetFieldName(s, f, &name);
2608: PetscSectionSetFieldName(sNew, f, name);
2609: PetscSectionGetFieldComponents(s, f, &numComp);
2610: PetscSectionSetFieldComponents(sNew, f, numComp);
2611: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2612: PetscSectionGetComponentName(s, f, c, &name);
2613: PetscSectionSetComponentName(sNew, f, c, name);
2614: }
2615: }
2616: ISGetLocalSize(permutation, &numPoints);
2617: ISGetIndices(permutation, &perm);
2618: PetscSectionGetChart(s, &pStart, &pEnd);
2619: PetscSectionSetChart(sNew, pStart, pEnd);
2620: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2621: for (p = pStart; p < pEnd; ++p) {
2622: PetscInt dof, cdof;
2624: PetscSectionGetDof(s, p, &dof);
2625: PetscSectionSetDof(sNew, perm[p], dof);
2626: PetscSectionGetConstraintDof(s, p, &cdof);
2627: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2628: for (f = 0; f < numFields; ++f) {
2629: PetscSectionGetFieldDof(s, p, f, &dof);
2630: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2631: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2632: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2633: }
2634: }
2635: PetscSectionSetUp(sNew);
2636: for (p = pStart; p < pEnd; ++p) {
2637: const PetscInt *cind;
2638: PetscInt cdof;
2640: PetscSectionGetConstraintDof(s, p, &cdof);
2641: if (cdof) {
2642: PetscSectionGetConstraintIndices(s, p, &cind);
2643: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2644: }
2645: for (f = 0; f < numFields; ++f) {
2646: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2647: if (cdof) {
2648: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2649: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2650: }
2651: }
2652: }
2653: ISRestoreIndices(permutation, &perm);
2654: *sectionNew = sNew;
2655: return(0);
2656: }
2658: /*@
2659: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2661: Collective on section
2663: Input Parameters:
2664: + section - The PetscSection
2665: . obj - A PetscObject which serves as the key for this index
2666: . clSection - Section giving the size of the closure of each point
2667: - clPoints - IS giving the points in each closure
2669: Note: We compress out closure points with no dofs in this section
2671: Level: advanced
2673: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2674: @*/
2675: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2676: {
2683: if (section->clObj != obj) {PetscSectionResetClosurePermutation(section);}
2684: section->clObj = obj;
2685: PetscObjectReference((PetscObject)clSection);
2686: PetscObjectReference((PetscObject)clPoints);
2687: PetscSectionDestroy(§ion->clSection);
2688: ISDestroy(§ion->clPoints);
2689: section->clSection = clSection;
2690: section->clPoints = clPoints;
2691: return(0);
2692: }
2694: /*@
2695: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2697: Collective on section
2699: Input Parameters:
2700: + section - The PetscSection
2701: - obj - A PetscObject which serves as the key for this index
2703: Output Parameters:
2704: + clSection - Section giving the size of the closure of each point
2705: - clPoints - IS giving the points in each closure
2707: Note: We compress out closure points with no dofs in this section
2709: Level: advanced
2711: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2712: @*/
2713: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2714: {
2716: if (section->clObj == obj) {
2717: if (clSection) *clSection = section->clSection;
2718: if (clPoints) *clPoints = section->clPoints;
2719: } else {
2720: if (clSection) *clSection = NULL;
2721: if (clPoints) *clPoints = NULL;
2722: }
2723: return(0);
2724: }
2726: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2727: {
2728: PetscInt i;
2729: khiter_t iter;
2730: int new_entry;
2731: PetscSectionClosurePermKey key = {depth, clSize};
2732: PetscSectionClosurePermVal *val;
2736: if (section->clObj != obj) {
2737: PetscSectionDestroy(§ion->clSection);
2738: ISDestroy(§ion->clPoints);
2739: }
2740: section->clObj = obj;
2741: if (!section->clHash) {PetscClPermCreate(§ion->clHash);}
2742: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
2743: val = &kh_val(section->clHash, iter);
2744: if (!new_entry) {
2745: PetscFree(val->perm);
2746: PetscFree(val->invPerm);
2747: }
2748: if (mode == PETSC_COPY_VALUES) {
2749: PetscMalloc1(clSize, &val->perm);
2750: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2751: PetscArraycpy(val->perm, clPerm, clSize);
2752: } else if (mode == PETSC_OWN_POINTER) {
2753: val->perm = clPerm;
2754: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2755: PetscMalloc1(clSize, &val->invPerm);
2756: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
2757: return(0);
2758: }
2760: /*@
2761: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2763: Not Collective
2765: Input Parameters:
2766: + section - The PetscSection
2767: . obj - A PetscObject which serves as the key for this index (usually a DM)
2768: . depth - Depth of points on which to apply the given permutation
2769: - perm - Permutation of the cell dof closure
2771: Note:
2772: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
2773: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
2774: each topology and degree.
2776: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
2777: exotic/enriched spaces on mixed topology meshes.
2779: Level: intermediate
2781: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2782: @*/
2783: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
2784: {
2785: const PetscInt *clPerm = NULL;
2786: PetscInt clSize = 0;
2787: PetscErrorCode ierr;
2790: if (perm) {
2791: ISGetLocalSize(perm, &clSize);
2792: ISGetIndices(perm, &clPerm);
2793: }
2794: PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2795: if (perm) {ISRestoreIndices(perm, &clPerm);}
2796: return(0);
2797: }
2799: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2800: {
2804: if (section->clObj == obj) {
2805: PetscSectionClosurePermKey k = {depth, size};
2806: PetscSectionClosurePermVal v;
2807: PetscClPermGet(section->clHash, k, &v);
2808: if (perm) *perm = v.perm;
2809: } else {
2810: if (perm) *perm = NULL;
2811: }
2812: return(0);
2813: }
2815: /*@
2816: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2818: Not collective
2820: Input Parameters:
2821: + section - The PetscSection
2822: . obj - A PetscObject which serves as the key for this index (usually a DM)
2823: . depth - Depth stratum on which to obtain closure permutation
2824: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2826: Output Parameter:
2827: . perm - The dof closure permutation
2829: Note:
2830: The user must destroy the IS that is returned.
2832: Level: intermediate
2834: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2835: @*/
2836: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2837: {
2838: const PetscInt *clPerm;
2839: PetscErrorCode ierr;
2842: PetscSectionGetClosurePermutation_Internal(section, obj, depth, clSize, &clPerm);
2843: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2844: return(0);
2845: }
2847: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
2848: {
2852: if (section->clObj == obj && section->clHash) {
2853: PetscSectionClosurePermKey k = {depth, size};
2854: PetscSectionClosurePermVal v;
2855: PetscClPermGet(section->clHash, k, &v);
2856: if (perm) *perm = v.invPerm;
2857: } else {
2858: if (perm) *perm = NULL;
2859: }
2860: return(0);
2861: }
2863: /*@
2864: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2866: Not collective
2868: Input Parameters:
2869: + section - The PetscSection
2870: . obj - A PetscObject which serves as the key for this index (usually a DM)
2871: . depth - Depth stratum on which to obtain closure permutation
2872: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
2874: Output Parameters:
2875: . perm - The dof closure permutation
2877: Note:
2878: The user must destroy the IS that is returned.
2880: Level: intermediate
2882: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2883: @*/
2884: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
2885: {
2886: const PetscInt *clPerm;
2887: PetscErrorCode ierr;
2890: PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm);
2891: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2892: return(0);
2893: }
2895: /*@
2896: PetscSectionGetField - Get the subsection associated with a single field
2898: Input Parameters:
2899: + s - The PetscSection
2900: - field - The field number
2902: Output Parameter:
2903: . subs - The subsection for the given field
2905: Level: intermediate
2907: .seealso: PetscSectionSetNumFields()
2908: @*/
2909: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2910: {
2914: 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);
2915: *subs = s->field[field];
2916: return(0);
2917: }
2919: PetscClassId PETSC_SECTION_SYM_CLASSID;
2920: PetscFunctionList PetscSectionSymList = NULL;
2922: /*@
2923: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
2925: Collective
2927: Input Parameter:
2928: . comm - the MPI communicator
2930: Output Parameter:
2931: . sym - pointer to the new set of symmetries
2933: Level: developer
2935: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2936: @*/
2937: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2938: {
2943: ISInitializePackage();
2944: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2945: return(0);
2946: }
2948: /*@C
2949: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
2951: Collective on PetscSectionSym
2953: Input Parameters:
2954: + sym - The section symmetry object
2955: - method - The name of the section symmetry type
2957: Level: developer
2959: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2960: @*/
2961: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2962: {
2963: PetscErrorCode (*r)(PetscSectionSym);
2964: PetscBool match;
2969: PetscObjectTypeCompare((PetscObject) sym, method, &match);
2970: if (match) return(0);
2972: PetscFunctionListFind(PetscSectionSymList,method,&r);
2973: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2974: if (sym->ops->destroy) {
2975: (*sym->ops->destroy)(sym);
2976: sym->ops->destroy = NULL;
2977: }
2978: (*r)(sym);
2979: PetscObjectChangeTypeName((PetscObject)sym,method);
2980: return(0);
2981: }
2983: /*@C
2984: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
2986: Not Collective
2988: Input Parameter:
2989: . sym - The section symmetry
2991: Output Parameter:
2992: . type - The index set type name
2994: Level: developer
2996: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2997: @*/
2998: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2999: {
3003: *type = ((PetscObject)sym)->type_name;
3004: return(0);
3005: }
3007: /*@C
3008: PetscSectionSymRegister - Adds a new section symmetry implementation
3010: Not Collective
3012: Input Parameters:
3013: + name - The name of a new user-defined creation routine
3014: - create_func - The creation routine itself
3016: Notes:
3017: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
3019: Level: developer
3021: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
3022: @*/
3023: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3024: {
3028: ISInitializePackage();
3029: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
3030: return(0);
3031: }
3033: /*@
3034: PetscSectionSymDestroy - Destroys a section symmetry.
3036: Collective on PetscSectionSym
3038: Input Parameters:
3039: . sym - the section symmetry
3041: Level: developer
3043: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
3044: @*/
3045: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3046: {
3047: SymWorkLink link,next;
3051: if (!*sym) return(0);
3053: if (--((PetscObject)(*sym))->refct > 0) {*sym = NULL; return(0);}
3054: if ((*sym)->ops->destroy) {
3055: (*(*sym)->ops->destroy)(*sym);
3056: }
3057: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3058: for (link=(*sym)->workin; link; link=next) {
3059: next = link->next;
3060: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
3061: PetscFree(link);
3062: }
3063: (*sym)->workin = NULL;
3064: PetscHeaderDestroy(sym);
3065: return(0);
3066: }
3068: /*@C
3069: PetscSectionSymView - Displays a section symmetry
3071: Collective on PetscSectionSym
3073: Input Parameters:
3074: + sym - the index set
3075: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
3077: Level: developer
3079: .seealso: PetscViewerASCIIOpen()
3080: @*/
3081: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3082: {
3087: if (!viewer) {
3088: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
3089: }
3092: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3093: if (sym->ops->view) {
3094: (*sym->ops->view)(sym,viewer);
3095: }
3096: return(0);
3097: }
3099: /*@
3100: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3102: Collective on PetscSection
3104: Input Parameters:
3105: + section - the section describing data layout
3106: - sym - the symmetry describing the affect of orientation on the access of the data
3108: Level: developer
3110: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3111: @*/
3112: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3113: {
3118: PetscSectionSymDestroy(&(section->sym));
3119: if (sym) {
3122: PetscObjectReference((PetscObject) sym);
3123: }
3124: section->sym = sym;
3125: return(0);
3126: }
3128: /*@
3129: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3131: Not collective
3133: Input Parameters:
3134: . section - the section describing data layout
3136: Output Parameters:
3137: . sym - the symmetry describing the affect of orientation on the access of the data
3139: Level: developer
3141: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3142: @*/
3143: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3144: {
3147: *sym = section->sym;
3148: return(0);
3149: }
3151: /*@
3152: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3154: Collective on PetscSection
3156: Input Parameters:
3157: + section - the section describing data layout
3158: . field - the field number
3159: - sym - the symmetry describing the affect of orientation on the access of the data
3161: Level: developer
3163: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3164: @*/
3165: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3166: {
3171: 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);
3172: PetscSectionSetSym(section->field[field],sym);
3173: return(0);
3174: }
3176: /*@
3177: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3179: Collective on PetscSection
3181: Input Parameters:
3182: + section - the section describing data layout
3183: - field - the field number
3185: Output Parameters:
3186: . sym - the symmetry describing the affect of orientation on the access of the data
3188: Level: developer
3190: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3191: @*/
3192: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3193: {
3196: 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);
3197: *sym = section->field[field]->sym;
3198: return(0);
3199: }
3201: /*@C
3202: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
3204: Not collective
3206: Input Parameters:
3207: + section - the section
3208: . numPoints - the number of points
3209: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3210: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3211: context, see DMPlexGetConeOrientation()).
3213: Output Parameters:
3214: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3215: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3216: identity).
3218: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3219: .vb
3220: const PetscInt **perms;
3221: const PetscScalar **rots;
3222: PetscInt lOffset;
3224: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3225: for (i = 0, lOffset = 0; i < numPoints; i++) {
3226: PetscInt point = points[2*i], dof, sOffset;
3227: const PetscInt *perm = perms ? perms[i] : NULL;
3228: const PetscScalar *rot = rots ? rots[i] : NULL;
3230: PetscSectionGetDof(section,point,&dof);
3231: PetscSectionGetOffset(section,point,&sOffset);
3233: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3234: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3235: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3236: lOffset += dof;
3237: }
3238: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3239: .ve
3241: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3242: .vb
3243: const PetscInt **perms;
3244: const PetscScalar **rots;
3245: PetscInt lOffset;
3247: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3248: for (i = 0, lOffset = 0; i < numPoints; i++) {
3249: PetscInt point = points[2*i], dof, sOffset;
3250: const PetscInt *perm = perms ? perms[i] : NULL;
3251: const PetscScalar *rot = rots ? rots[i] : NULL;
3253: PetscSectionGetDof(section,point,&dof);
3254: PetscSectionGetOffset(section,point,&sOff);
3256: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3257: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3258: offset += dof;
3259: }
3260: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3261: .ve
3263: Level: developer
3265: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3266: @*/
3267: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3268: {
3269: PetscSectionSym sym;
3270: PetscErrorCode ierr;
3275: if (perms) *perms = NULL;
3276: if (rots) *rots = NULL;
3277: sym = section->sym;
3278: if (sym && (perms || rots)) {
3279: SymWorkLink link;
3281: if (sym->workin) {
3282: link = sym->workin;
3283: sym->workin = sym->workin->next;
3284: } else {
3285: PetscNewLog(sym,&link);
3286: }
3287: if (numPoints > link->numPoints) {
3288: PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3289: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3290: link->numPoints = numPoints;
3291: }
3292: link->next = sym->workout;
3293: sym->workout = link;
3294: PetscArrayzero((PetscInt**)link->perms,numPoints);
3295: PetscArrayzero((PetscInt**)link->rots,numPoints);
3296: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3297: if (perms) *perms = link->perms;
3298: if (rots) *rots = link->rots;
3299: }
3300: return(0);
3301: }
3303: /*@C
3304: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
3306: Not collective
3308: Input Parameters:
3309: + section - the section
3310: . numPoints - the number of points
3311: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3312: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3313: context, see DMPlexGetConeOrientation()).
3315: Output Parameters:
3316: + perms - The permutations for the given orientations: set to NULL at conclusion
3317: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3319: Level: developer
3321: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3322: @*/
3323: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3324: {
3325: PetscSectionSym sym;
3329: sym = section->sym;
3330: if (sym && (perms || rots)) {
3331: SymWorkLink *p,link;
3333: for (p=&sym->workout; (link=*p); p=&link->next) {
3334: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3335: *p = link->next;
3336: link->next = sym->workin;
3337: sym->workin = link;
3338: if (perms) *perms = NULL;
3339: if (rots) *rots = NULL;
3340: return(0);
3341: }
3342: }
3343: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3344: }
3345: return(0);
3346: }
3348: /*@C
3349: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
3351: Not collective
3353: Input Parameters:
3354: + section - the section
3355: . field - the field of the section
3356: . numPoints - the number of points
3357: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3358: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3359: context, see DMPlexGetConeOrientation()).
3361: Output Parameters:
3362: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3363: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3364: identity).
3366: Level: developer
3368: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3369: @*/
3370: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3371: {
3376: 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);
3377: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3378: return(0);
3379: }
3381: /*@C
3382: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
3384: Not collective
3386: Input Parameters:
3387: + section - the section
3388: . field - the field number
3389: . numPoints - the number of points
3390: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3391: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3392: context, see DMPlexGetConeOrientation()).
3394: Output Parameters:
3395: + perms - The permutations for the given orientations: set to NULL at conclusion
3396: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3398: Level: developer
3400: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3401: @*/
3402: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3403: {
3408: 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);
3409: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3410: return(0);
3411: }
3413: /*@
3414: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3416: Not collective
3418: Input Parameter:
3419: . s - the global PetscSection
3421: Output Parameters:
3422: . flg - the flag
3424: Level: developer
3426: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3427: @*/
3428: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3429: {
3432: *flg = s->useFieldOff;
3433: return(0);
3434: }
3436: /*@
3437: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3439: Not collective
3441: Input Parameters:
3442: + s - the global PetscSection
3443: - flg - the flag
3445: Level: developer
3447: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3448: @*/
3449: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3450: {
3453: s->useFieldOff = flg;
3454: return(0);
3455: }
3457: #define PetscSectionExpandPoints_Loop(TYPE) \
3458: { \
3459: PetscInt i, n, o0, o1, size; \
3460: TYPE *a0 = (TYPE*)origArray, *a1; \
3461: PetscSectionGetStorageSize(s, &size); \
3462: PetscMalloc1(size, &a1); \
3463: for (i=0; i<npoints; i++) { \
3464: PetscSectionGetOffset(origSection, points_[i], &o0); \
3465: PetscSectionGetOffset(s, i, &o1); \
3466: PetscSectionGetDof(s, i, &n); \
3467: PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3468: } \
3469: *newArray = (void*)a1; \
3470: }
3472: /*@
3473: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3475: Not collective
3477: Input Parameters:
3478: + origSection - the PetscSection describing the layout of the array
3479: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3480: . origArray - the array; its size must be equal to the storage size of origSection
3481: - points - IS with points to extract; its indices must lie in the chart of origSection
3483: Output Parameters:
3484: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3485: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3487: Level: developer
3489: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3490: @*/
3491: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3492: {
3493: PetscSection s;
3494: const PetscInt *points_;
3495: PetscInt i, n, npoints, pStart, pEnd;
3496: PetscMPIInt unitsize;
3497: PetscErrorCode ierr;
3505: MPI_Type_size(dataType, &unitsize);
3506: ISGetLocalSize(points, &npoints);
3507: ISGetIndices(points, &points_);
3508: PetscSectionGetChart(origSection, &pStart, &pEnd);
3509: PetscSectionCreate(PETSC_COMM_SELF, &s);
3510: PetscSectionSetChart(s, 0, npoints);
3511: for (i=0; i<npoints; i++) {
3512: 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);
3513: PetscSectionGetDof(origSection, points_[i], &n);
3514: PetscSectionSetDof(s, i, n);
3515: }
3516: PetscSectionSetUp(s);
3517: if (newArray) {
3518: if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);}
3519: else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);}
3520: else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);}
3521: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3522: }
3523: if (newSection) {
3524: *newSection = s;
3525: } else {
3526: PetscSectionDestroy(&s);
3527: }
3528: ISRestoreIndices(points, &points_);
3529: return(0);
3530: }