Actual source code: vsectionis.c
petsc-3.11.4 2019-09-28
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/isimpl.h>
6: #include <petscsf.h>
7: #include <petscviewer.h>
9: PetscClassId PETSC_SECTION_CLASSID;
11: /*@
12: PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
14: Collective on MPI_Comm
16: Input Parameters:
17: + comm - the MPI communicator
18: - s - pointer to the section
20: Level: developer
22: Notes:
23: Typical calling sequence
24: $ PetscSectionCreate(MPI_Comm,PetscSection *);
25: $ PetscSectionSetNumFields(PetscSection, numFields);
26: $ PetscSectionSetChart(PetscSection,low,high);
27: $ PetscSectionSetDof(PetscSection,point,numdof);
28: $ PetscSectionSetUp(PetscSection);
29: $ PetscSectionGetOffset(PetscSection,point,PetscInt *);
30: $ PetscSectionDestroy(PetscSection);
32: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
33: recommended they not be used in user codes unless you really gain something in their use.
35: .seealso: PetscSection, PetscSectionDestroy()
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
43: ISInitializePackage();
45: PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
47: (*s)->pStart = -1;
48: (*s)->pEnd = -1;
49: (*s)->perm = NULL;
50: (*s)->pointMajor = 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)->clObj = NULL;
62: (*s)->clSection = NULL;
63: (*s)->clPoints = NULL;
64: (*s)->clSize = 0;
65: (*s)->clPerm = NULL;
66: (*s)->clInvPerm = NULL;
67: return(0);
68: }
70: /*@
71: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
73: Collective on MPI_Comm
75: Input Parameter:
76: . section - the PetscSection
78: Output Parameter:
79: . newSection - the copy
81: Level: developer
83: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
84: @*/
85: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
86: {
87: PetscSectionSym sym;
88: IS perm;
89: PetscInt numFields, f, pStart, pEnd, p;
90: PetscErrorCode ierr;
95: PetscSectionGetNumFields(section, &numFields);
96: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
97: for (f = 0; f < numFields; ++f) {
98: const char *name = NULL;
99: PetscInt numComp = 0;
101: PetscSectionGetFieldName(section, f, &name);
102: PetscSectionSetFieldName(newSection, f, name);
103: PetscSectionGetFieldComponents(section, f, &numComp);
104: PetscSectionSetFieldComponents(newSection, f, numComp);
105: PetscSectionGetFieldSym(section, f, &sym);
106: PetscSectionSetFieldSym(newSection, f, sym);
107: }
108: PetscSectionGetChart(section, &pStart, &pEnd);
109: PetscSectionSetChart(newSection, pStart, pEnd);
110: PetscSectionGetPermutation(section, &perm);
111: PetscSectionSetPermutation(newSection, perm);
112: PetscSectionGetSym(section, &sym);
113: PetscSectionSetSym(newSection, sym);
114: for (p = pStart; p < pEnd; ++p) {
115: PetscInt dof, cdof, fcdof = 0;
117: PetscSectionGetDof(section, p, &dof);
118: PetscSectionSetDof(newSection, p, dof);
119: PetscSectionGetConstraintDof(section, p, &cdof);
120: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
121: for (f = 0; f < numFields; ++f) {
122: PetscSectionGetFieldDof(section, p, f, &dof);
123: PetscSectionSetFieldDof(newSection, p, f, dof);
124: if (cdof) {
125: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
126: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
127: }
128: }
129: }
130: PetscSectionSetUp(newSection);
131: for (p = pStart; p < pEnd; ++p) {
132: PetscInt off, cdof, fcdof = 0;
133: const PetscInt *cInd;
135: /* Must set offsets in case they do not agree with the prefix sums */
136: PetscSectionGetOffset(section, p, &off);
137: PetscSectionSetOffset(newSection, p, off);
138: PetscSectionGetConstraintDof(section, p, &cdof);
139: if (cdof) {
140: PetscSectionGetConstraintIndices(section, p, &cInd);
141: PetscSectionSetConstraintIndices(newSection, p, cInd);
142: for (f = 0; f < numFields; ++f) {
143: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
144: if (fcdof) {
145: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
146: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
147: }
148: }
149: }
150: }
151: return(0);
152: }
154: /*@
155: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
157: Collective on MPI_Comm
159: Input Parameter:
160: . section - the PetscSection
162: Output Parameter:
163: . newSection - the copy
165: Level: developer
167: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
168: @*/
169: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
170: {
176: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
177: PetscSectionCopy(section, *newSection);
178: return(0);
179: }
181: /*@
182: PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database
184: Collective on PetscSection
186: Input Parameter:
187: . section - the PetscSection
189: Options Database:
190: . -petscsection_point_major the dof order
192: Level: intermediate
194: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
195: @*/
196: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
197: {
202: PetscObjectOptionsBegin((PetscObject) s);
203: PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
204: /* process any options handlers added with PetscObjectAddOptionsHandler() */
205: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
206: PetscOptionsEnd();
207: PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
208: return(0);
209: }
211: /*@
212: PetscSectionCompare - Compares two sections
214: Collective on PetscSection
216: Input Parameterss:
217: + s1 - the first PetscSection
218: - s2 - the second PetscSection
220: Output Parameter:
221: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
223: Level: developer
225: Notes:
226: Field names are disregarded.
228: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
229: @*/
230: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
231: {
232: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
233: const PetscInt *idx1, *idx2;
234: IS perm1, perm2;
235: PetscBool flg;
236: PetscMPIInt mflg;
243: flg = PETSC_FALSE;
245: MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
246: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
247: *congruent = PETSC_FALSE;
248: return(0);
249: }
251: PetscSectionGetChart(s1, &pStart, &pEnd);
252: PetscSectionGetChart(s2, &n1, &n2);
253: if (pStart != n1 || pEnd != n2) goto not_congruent;
255: PetscSectionGetPermutation(s1, &perm1);
256: PetscSectionGetPermutation(s2, &perm2);
257: if (perm1 && perm2) {
258: ISEqual(perm1, perm2, congruent);
259: if (!(*congruent)) goto not_congruent;
260: } else if (perm1 != perm2) goto not_congruent;
262: for (p = pStart; p < pEnd; ++p) {
263: PetscSectionGetOffset(s1, p, &n1);
264: PetscSectionGetOffset(s2, p, &n2);
265: if (n1 != n2) goto not_congruent;
267: PetscSectionGetDof(s1, p, &n1);
268: PetscSectionGetDof(s2, p, &n2);
269: if (n1 != n2) goto not_congruent;
271: PetscSectionGetConstraintDof(s1, p, &ncdof);
272: PetscSectionGetConstraintDof(s2, p, &n2);
273: if (ncdof != n2) goto not_congruent;
275: PetscSectionGetConstraintIndices(s1, p, &idx1);
276: PetscSectionGetConstraintIndices(s2, p, &idx2);
277: PetscMemcmp(idx1, idx2, ncdof*sizeof(PetscInt), congruent);
278: if (!(*congruent)) goto not_congruent;
279: }
281: PetscSectionGetNumFields(s1, &nfields);
282: PetscSectionGetNumFields(s2, &n2);
283: if (nfields != n2) goto not_congruent;
285: for (f = 0; f < nfields; ++f) {
286: PetscSectionGetFieldComponents(s1, f, &n1);
287: PetscSectionGetFieldComponents(s2, f, &n2);
288: if (n1 != n2) goto not_congruent;
290: for (p = pStart; p < pEnd; ++p) {
291: PetscSectionGetFieldOffset(s1, p, f, &n1);
292: PetscSectionGetFieldOffset(s2, p, f, &n2);
293: if (n1 != n2) goto not_congruent;
295: PetscSectionGetFieldDof(s1, p, f, &n1);
296: PetscSectionGetFieldDof(s2, p, f, &n2);
297: if (n1 != n2) goto not_congruent;
299: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
300: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
301: if (nfcdof != n2) goto not_congruent;
303: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
304: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
305: PetscMemcmp(idx1, idx2, nfcdof*sizeof(PetscInt), congruent);
306: if (!(*congruent)) goto not_congruent;
307: }
308: }
310: flg = PETSC_TRUE;
311: not_congruent:
312: MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
313: return(0);
314: }
316: /*@
317: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
319: Not collective
321: Input Parameter:
322: . s - the PetscSection
324: Output Parameter:
325: . numFields - the number of fields defined, or 0 if none were defined
327: Level: intermediate
329: .seealso: PetscSectionSetNumFields()
330: @*/
331: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
332: {
336: *numFields = s->numFields;
337: return(0);
338: }
340: /*@
341: PetscSectionSetNumFields - Sets the number of fields.
343: Not collective
345: Input Parameters:
346: + s - the PetscSection
347: - numFields - the number of fields
349: Level: intermediate
351: .seealso: PetscSectionGetNumFields()
352: @*/
353: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
354: {
355: PetscInt f;
360: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %D must be positive", numFields);
361: PetscSectionReset(s);
363: s->numFields = numFields;
364: PetscMalloc1(s->numFields, &s->numFieldComponents);
365: PetscMalloc1(s->numFields, &s->fieldNames);
366: PetscMalloc1(s->numFields, &s->field);
367: for (f = 0; f < s->numFields; ++f) {
368: char name[64];
370: s->numFieldComponents[f] = 1;
372: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
373: PetscSNPrintf(name, 64, "Field_%D", f);
374: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
375: }
376: return(0);
377: }
379: /*@C
380: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
382: Not Collective
384: Input Parameters:
385: + s - the PetscSection
386: - field - the field number
388: Output Parameter:
389: . fieldName - the field name
391: Level: developer
393: .seealso: PetscSectionSetFieldName()
394: @*/
395: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
396: {
400: 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);
401: *fieldName = s->fieldNames[field];
402: return(0);
403: }
405: /*@C
406: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
408: Not Collective
410: Input Parameters:
411: + s - the PetscSection
412: . field - the field number
413: - fieldName - the field name
415: Level: developer
417: .seealso: PetscSectionGetFieldName()
418: @*/
419: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
420: {
426: 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);
427: PetscFree(s->fieldNames[field]);
428: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
429: return(0);
430: }
432: /*@
433: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
435: Not collective
437: Input Parameters:
438: + s - the PetscSection
439: - field - the field number
441: Output Parameter:
442: . numComp - the number of field components
444: Level: intermediate
446: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
447: @*/
448: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
449: {
453: 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);
454: *numComp = s->numFieldComponents[field];
455: return(0);
456: }
458: /*@
459: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
461: Not collective
463: Input Parameters:
464: + s - the PetscSection
465: . field - the field number
466: - numComp - the number of field components
468: Level: intermediate
470: .seealso: PetscSectionGetNumFieldComponents(), PetscSectionGetNumFields()
471: @*/
472: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
473: {
476: 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);
477: s->numFieldComponents[field] = numComp;
478: return(0);
479: }
481: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
482: {
486: if (!s->bc) {
487: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
488: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
489: }
490: return(0);
491: }
493: /*@
494: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
496: Not collective
498: Input Parameter:
499: . s - the PetscSection
501: Output Parameters:
502: + pStart - the first point
503: - pEnd - one past the last point
505: Level: intermediate
507: .seealso: PetscSectionSetChart(), PetscSectionCreate()
508: @*/
509: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
510: {
513: if (pStart) *pStart = s->pStart;
514: if (pEnd) *pEnd = s->pEnd;
515: return(0);
516: }
518: /*@
519: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
521: Not collective
523: Input Parameters:
524: + s - the PetscSection
525: . pStart - the first point
526: - pEnd - one past the last point
528: Level: intermediate
530: .seealso: PetscSectionGetChart(), PetscSectionCreate()
531: @*/
532: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
533: {
534: PetscInt f;
539: /* Cannot Reset() because it destroys field information */
540: s->setup = PETSC_FALSE;
541: PetscSectionDestroy(&s->bc);
542: PetscFree(s->bcIndices);
543: PetscFree2(s->atlasDof, s->atlasOff);
545: s->pStart = pStart;
546: s->pEnd = pEnd;
547: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
548: PetscMemzero(s->atlasDof, (pEnd - pStart)*sizeof(PetscInt));
549: for (f = 0; f < s->numFields; ++f) {
550: PetscSectionSetChart(s->field[f], pStart, pEnd);
551: }
552: return(0);
553: }
555: /*@
556: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
558: Not collective
560: Input Parameter:
561: . s - the PetscSection
563: Output Parameters:
564: . perm - The permutation as an IS
566: Level: intermediate
568: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
569: @*/
570: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
571: {
575: return(0);
576: }
578: /*@
579: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
581: Not collective
583: Input Parameters:
584: + s - the PetscSection
585: - perm - the permutation of points
587: Level: intermediate
589: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
590: @*/
591: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
592: {
598: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
599: if (s->perm != perm) {
600: ISDestroy(&s->perm);
601: if (perm) {
602: s->perm = perm;
603: PetscObjectReference((PetscObject) s->perm);
604: }
605: }
606: return(0);
607: }
609: /*@
610: PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major
612: Not collective
614: Input Parameter:
615: . s - the PetscSection
617: Output Parameter:
618: . pm - the flag for point major ordering
620: Level: intermediate
622: .seealso: PetscSectionSetPointMajor()
623: @*/
624: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
625: {
629: *pm = s->pointMajor;
630: return(0);
631: }
633: /*@
634: PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major
636: Not collective
638: Input Parameters:
639: + s - the PetscSection
640: - pm - the flag for point major ordering
642: Not collective
644: Level: intermediate
646: .seealso: PetscSectionGetPointMajor()
647: @*/
648: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
649: {
652: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
653: s->pointMajor = pm;
654: return(0);
655: }
657: /*@
658: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
660: Not collective
662: Input Parameters:
663: + s - the PetscSection
664: - point - the point
666: Output Parameter:
667: . numDof - the number of dof
669: Level: intermediate
671: .seealso: PetscSectionSetDof(), PetscSectionCreate()
672: @*/
673: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
674: {
678: #if defined(PETSC_USE_DEBUG)
679: 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);
680: #endif
681: *numDof = s->atlasDof[point - s->pStart];
682: return(0);
683: }
685: /*@
686: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
688: Not collective
690: Input Parameters:
691: + s - the PetscSection
692: . point - the point
693: - numDof - the number of dof
695: Level: intermediate
697: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
698: @*/
699: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
700: {
703: 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);
704: s->atlasDof[point - s->pStart] = numDof;
705: return(0);
706: }
708: /*@
709: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
711: Not collective
713: Input Parameters:
714: + s - the PetscSection
715: . point - the point
716: - numDof - the number of additional dof
718: Level: intermediate
720: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
721: @*/
722: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
723: {
726: 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);
727: s->atlasDof[point - s->pStart] += numDof;
728: return(0);
729: }
731: /*@
732: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
734: Not collective
736: Input Parameters:
737: + s - the PetscSection
738: . point - the point
739: - field - the field
741: Output Parameter:
742: . numDof - the number of dof
744: Level: intermediate
746: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
747: @*/
748: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
749: {
754: 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);
755: PetscSectionGetDof(s->field[field], point, numDof);
756: return(0);
757: }
759: /*@
760: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
762: Not collective
764: Input Parameters:
765: + s - the PetscSection
766: . point - the point
767: . field - the field
768: - numDof - the number of dof
770: Level: intermediate
772: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
773: @*/
774: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
775: {
780: 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);
781: PetscSectionSetDof(s->field[field], point, numDof);
782: return(0);
783: }
785: /*@
786: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
788: Not collective
790: Input Parameters:
791: + s - the PetscSection
792: . point - the point
793: . field - the field
794: - numDof - the number of dof
796: Level: intermediate
798: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
799: @*/
800: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
801: {
806: 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);
807: PetscSectionAddDof(s->field[field], point, numDof);
808: return(0);
809: }
811: /*@
812: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
814: Not collective
816: Input Parameters:
817: + s - the PetscSection
818: - point - the point
820: Output Parameter:
821: . numDof - the number of dof which are fixed by constraints
823: Level: intermediate
825: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
826: @*/
827: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
828: {
834: if (s->bc) {
835: PetscSectionGetDof(s->bc, point, numDof);
836: } else *numDof = 0;
837: return(0);
838: }
840: /*@
841: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
843: Not collective
845: Input Parameters:
846: + s - the PetscSection
847: . point - the point
848: - numDof - the number of dof which are fixed by constraints
850: Level: intermediate
852: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
853: @*/
854: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
855: {
860: if (numDof) {
861: PetscSectionCheckConstraints_Static(s);
862: PetscSectionSetDof(s->bc, point, numDof);
863: }
864: return(0);
865: }
867: /*@
868: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
870: Not collective
872: Input Parameters:
873: + s - the PetscSection
874: . point - the point
875: - numDof - the number of additional dof which are fixed by constraints
877: Level: intermediate
879: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
880: @*/
881: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
882: {
887: if (numDof) {
888: PetscSectionCheckConstraints_Static(s);
889: PetscSectionAddDof(s->bc, point, numDof);
890: }
891: return(0);
892: }
894: /*@
895: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
897: Not collective
899: Input Parameters:
900: + s - the PetscSection
901: . point - the point
902: - field - the field
904: Output Parameter:
905: . numDof - the number of dof which are fixed by constraints
907: Level: intermediate
909: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
910: @*/
911: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
912: {
918: 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);
919: PetscSectionGetConstraintDof(s->field[field], point, numDof);
920: return(0);
921: }
923: /*@
924: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
926: Not collective
928: Input Parameters:
929: + s - the PetscSection
930: . point - the point
931: . field - the field
932: - numDof - the number of dof which are fixed by constraints
934: Level: intermediate
936: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
937: @*/
938: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
939: {
944: 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);
945: PetscSectionSetConstraintDof(s->field[field], point, numDof);
946: return(0);
947: }
949: /*@
950: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
952: Not collective
954: Input Parameters:
955: + s - the PetscSection
956: . point - the point
957: . field - the field
958: - numDof - the number of additional dof which are fixed by constraints
960: Level: intermediate
962: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
963: @*/
964: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
965: {
970: 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);
971: PetscSectionAddConstraintDof(s->field[field], point, numDof);
972: return(0);
973: }
975: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
976: {
981: if (s->bc) {
982: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
984: PetscSectionSetUp(s->bc);
985: PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices);
986: }
987: return(0);
988: }
990: /*@
991: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
993: Not collective
995: Input Parameter:
996: . s - the PetscSection
998: Level: intermediate
1000: .seealso: PetscSectionCreate()
1001: @*/
1002: PetscErrorCode PetscSectionSetUp(PetscSection s)
1003: {
1004: const PetscInt *pind = NULL;
1005: PetscInt offset = 0, foff, p, f;
1006: PetscErrorCode ierr;
1010: if (s->setup) return(0);
1011: s->setup = PETSC_TRUE;
1012: /* Set offsets and field offsets for all points */
1013: /* Assume that all fields have the same chart */
1014: if (s->perm) {ISGetIndices(s->perm, &pind);}
1015: if (s->pointMajor) {
1016: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1017: const PetscInt q = pind ? pind[p] : p;
1019: /* Set point offset */
1020: s->atlasOff[q] = offset;
1021: offset += s->atlasDof[q];
1022: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
1023: /* Set field offset */
1024: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1025: PetscSection sf = s->field[f];
1027: sf->atlasOff[q] = foff;
1028: foff += sf->atlasDof[q];
1029: }
1030: }
1031: } else {
1032: /* Set field offsets for all points */
1033: for (f = 0; f < s->numFields; ++f) {
1034: PetscSection sf = s->field[f];
1036: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1037: const PetscInt q = pind ? pind[p] : p;
1039: sf->atlasOff[q] = offset;
1040: offset += sf->atlasDof[q];
1041: }
1042: }
1043: /* Disable point offsets since these are unused */
1044: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1045: s->atlasOff[p] = -1;
1046: s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1047: }
1048: }
1049: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1050: /* Setup BC sections */
1051: PetscSectionSetUpBC(s);
1052: for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1053: return(0);
1054: }
1056: /*@
1057: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1059: Not collective
1061: Input Parameters:
1062: . s - the PetscSection
1064: Output Parameter:
1065: . maxDof - the maximum dof
1067: Level: intermediate
1069: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1070: @*/
1071: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1072: {
1076: *maxDof = s->maxDof;
1077: return(0);
1078: }
1080: /*@
1081: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
1083: Not collective
1085: Input Parameters:
1086: + s - the PetscSection
1087: - size - the allocated size
1089: Output Parameter:
1090: . size - the size of an array which can hold all the dofs
1092: Level: intermediate
1094: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1095: @*/
1096: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1097: {
1098: PetscInt p, n = 0;
1103: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1104: *size = n;
1105: return(0);
1106: }
1108: /*@
1109: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1111: Not collective
1113: Input Parameters:
1114: + s - the PetscSection
1115: - point - the point
1117: Output Parameter:
1118: . size - the size of an array which can hold all unconstrained dofs
1120: Level: intermediate
1122: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1123: @*/
1124: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1125: {
1126: PetscInt p, n = 0;
1131: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1132: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1133: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1134: }
1135: *size = n;
1136: return(0);
1137: }
1139: /*@
1140: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1141: the local section and an SF describing the section point overlap.
1143: Input Parameters:
1144: + s - The PetscSection for the local field layout
1145: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1146: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1147: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1149: Output Parameter:
1150: . gsection - The PetscSection for the global field layout
1152: Note: This gives negative sizes and offsets to points not owned by this process
1154: Level: developer
1156: .seealso: PetscSectionCreate()
1157: @*/
1158: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1159: {
1160: PetscSection gs;
1161: const PetscInt *pind = NULL;
1162: PetscInt *recv = NULL, *neg = NULL;
1163: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal;
1164: PetscErrorCode ierr;
1172: PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1173: PetscSectionGetChart(s, &pStart, &pEnd);
1174: PetscSectionSetChart(gs, pStart, pEnd);
1175: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1176: nlocal = nroots; /* The local/leaf space matches global/root space */
1177: /* Must allocate for all points visible to SF, which may be more than this section */
1178: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1179: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1180: PetscMalloc2(nroots,&neg,nlocal,&recv);
1181: PetscMemzero(neg,nroots*sizeof(PetscInt));
1182: }
1183: /* Mark all local points with negative dof */
1184: for (p = pStart; p < pEnd; ++p) {
1185: PetscSectionGetDof(s, p, &dof);
1186: PetscSectionSetDof(gs, p, dof);
1187: PetscSectionGetConstraintDof(s, p, &cdof);
1188: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1189: if (neg) neg[p] = -(dof+1);
1190: }
1191: PetscSectionSetUpBC(gs);
1192: if (gs->bcIndices) {PetscMemcpy(gs->bcIndices, s->bcIndices, sizeof(PetscInt) * (gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]));}
1193: if (nroots >= 0) {
1194: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1195: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1196: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1197: for (p = pStart; p < pEnd; ++p) {
1198: if (recv[p] < 0) {
1199: gs->atlasDof[p-pStart] = recv[p];
1200: PetscSectionGetDof(s, p, &dof);
1201: 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);
1202: }
1203: }
1204: }
1205: /* Calculate new sizes, get process offset, and calculate point offsets */
1206: if (s->perm) {ISGetIndices(s->perm, &pind);}
1207: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1208: const PetscInt q = pind ? pind[p] : p;
1210: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1211: gs->atlasOff[q] = off;
1212: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1213: }
1214: if (!localOffsets) {
1215: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1216: globalOff -= off;
1217: }
1218: for (p = pStart, off = 0; p < pEnd; ++p) {
1219: gs->atlasOff[p-pStart] += globalOff;
1220: if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1221: }
1222: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1223: /* Put in negative offsets for ghost points */
1224: if (nroots >= 0) {
1225: PetscMemzero(recv,nlocal*sizeof(PetscInt));
1226: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1227: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1228: for (p = pStart; p < pEnd; ++p) {
1229: if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1230: }
1231: }
1232: PetscFree2(neg,recv);
1233: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1234: *gsection = gs;
1235: return(0);
1236: }
1238: /*@
1239: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1240: the local section and an SF describing the section point overlap.
1242: Input Parameters:
1243: + s - The PetscSection for the local field layout
1244: . sf - The SF describing parallel layout of the section points
1245: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1246: . numExcludes - The number of exclusion ranges
1247: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1249: Output Parameter:
1250: . gsection - The PetscSection for the global field layout
1252: Note: This gives negative sizes and offsets to points not owned by this process
1254: Level: developer
1256: .seealso: PetscSectionCreate()
1257: @*/
1258: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1259: {
1260: const PetscInt *pind = NULL;
1261: PetscInt *neg = NULL, *tmpOff = NULL;
1262: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1263: PetscErrorCode ierr;
1269: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1270: PetscSectionGetChart(s, &pStart, &pEnd);
1271: PetscSectionSetChart(*gsection, pStart, pEnd);
1272: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1273: if (nroots >= 0) {
1274: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1275: PetscCalloc1(nroots, &neg);
1276: if (nroots > pEnd-pStart) {
1277: PetscCalloc1(nroots, &tmpOff);
1278: } else {
1279: tmpOff = &(*gsection)->atlasDof[-pStart];
1280: }
1281: }
1282: /* Mark ghost points with negative dof */
1283: for (p = pStart; p < pEnd; ++p) {
1284: for (e = 0; e < numExcludes; ++e) {
1285: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1286: PetscSectionSetDof(*gsection, p, 0);
1287: break;
1288: }
1289: }
1290: if (e < numExcludes) continue;
1291: PetscSectionGetDof(s, p, &dof);
1292: PetscSectionSetDof(*gsection, p, dof);
1293: PetscSectionGetConstraintDof(s, p, &cdof);
1294: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1295: if (neg) neg[p] = -(dof+1);
1296: }
1297: PetscSectionSetUpBC(*gsection);
1298: if (nroots >= 0) {
1299: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1300: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1301: if (nroots > pEnd - pStart) {
1302: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1303: }
1304: }
1305: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1306: if (s->perm) {ISGetIndices(s->perm, &pind);}
1307: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1308: const PetscInt q = pind ? pind[p] : p;
1310: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1311: (*gsection)->atlasOff[q] = off;
1312: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1313: }
1314: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1315: globalOff -= off;
1316: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1317: (*gsection)->atlasOff[p] += globalOff;
1318: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1319: }
1320: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1321: /* Put in negative offsets for ghost points */
1322: if (nroots >= 0) {
1323: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1324: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1325: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1326: if (nroots > pEnd - pStart) {
1327: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1328: }
1329: }
1330: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1331: PetscFree(neg);
1332: return(0);
1333: }
1335: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1336: {
1337: PetscInt pStart, pEnd, p, localSize = 0;
1341: PetscSectionGetChart(s, &pStart, &pEnd);
1342: for (p = pStart; p < pEnd; ++p) {
1343: PetscInt dof;
1345: PetscSectionGetDof(s, p, &dof);
1346: if (dof > 0) ++localSize;
1347: }
1348: PetscLayoutCreate(comm, layout);
1349: PetscLayoutSetLocalSize(*layout, localSize);
1350: PetscLayoutSetBlockSize(*layout, 1);
1351: PetscLayoutSetUp(*layout);
1352: return(0);
1353: }
1355: /*@
1356: PetscSectionGetValueLayout - Get the PetscLayout associated with a section, usually the default global section.
1358: Input Parameters:
1359: + comm - The MPI_Comm
1360: - s - The PetscSection
1362: Output Parameter:
1363: . layout - The layout for the section
1365: Level: developer
1367: .seealso: PetscSectionCreate()
1368: @*/
1369: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1370: {
1371: PetscInt pStart, pEnd, p, localSize = 0;
1377: PetscSectionGetChart(s, &pStart, &pEnd);
1378: for (p = pStart; p < pEnd; ++p) {
1379: PetscInt dof,cdof;
1381: PetscSectionGetDof(s, p, &dof);
1382: PetscSectionGetConstraintDof(s, p, &cdof);
1383: if (dof-cdof > 0) localSize += dof-cdof;
1384: }
1385: PetscLayoutCreate(comm, layout);
1386: PetscLayoutSetLocalSize(*layout, localSize);
1387: PetscLayoutSetBlockSize(*layout, 1);
1388: PetscLayoutSetUp(*layout);
1389: return(0);
1390: }
1392: /*@
1393: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1395: Not collective
1397: Input Parameters:
1398: + s - the PetscSection
1399: - point - the point
1401: Output Parameter:
1402: . offset - the offset
1404: Level: intermediate
1406: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1407: @*/
1408: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1409: {
1413: #if defined(PETSC_USE_DEBUG)
1414: 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);
1415: #endif
1416: *offset = s->atlasOff[point - s->pStart];
1417: return(0);
1418: }
1420: /*@
1421: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1423: Not collective
1425: Input Parameters:
1426: + s - the PetscSection
1427: . point - the point
1428: - offset - the offset
1430: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1432: Level: intermediate
1434: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1435: @*/
1436: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1437: {
1440: 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);
1441: s->atlasOff[point - s->pStart] = offset;
1442: return(0);
1443: }
1445: /*@
1446: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1448: Not collective
1450: Input Parameters:
1451: + s - the PetscSection
1452: . point - the point
1453: - field - the field
1455: Output Parameter:
1456: . offset - the offset
1458: Level: intermediate
1460: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1461: @*/
1462: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1463: {
1469: 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);
1470: PetscSectionGetOffset(s->field[field], point, offset);
1471: return(0);
1472: }
1474: /*@
1475: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1477: Not collective
1479: Input Parameters:
1480: + s - the PetscSection
1481: . point - the point
1482: . field - the field
1483: - offset - the offset
1485: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1487: Level: intermediate
1489: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1490: @*/
1491: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1492: {
1497: 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);
1498: PetscSectionSetOffset(s->field[field], point, offset);
1499: return(0);
1500: }
1502: /* This gives the offset on a point of the field, ignoring constraints */
1503: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1504: {
1505: PetscInt off, foff;
1511: 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);
1512: PetscSectionGetOffset(s, point, &off);
1513: PetscSectionGetOffset(s->field[field], point, &foff);
1514: *offset = foff - off;
1515: return(0);
1516: }
1518: /*@
1519: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1521: Not collective
1523: Input Parameter:
1524: . s - the PetscSection
1526: Output Parameters:
1527: + start - the minimum offset
1528: - end - one more than the maximum offset
1530: Level: intermediate
1532: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1533: @*/
1534: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1535: {
1536: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1541: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1542: PetscSectionGetChart(s, &pStart, &pEnd);
1543: for (p = 0; p < pEnd-pStart; ++p) {
1544: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1546: if (off >= 0) {
1547: os = PetscMin(os, off);
1548: oe = PetscMax(oe, off+dof);
1549: }
1550: }
1551: if (start) *start = os;
1552: if (end) *end = oe;
1553: return(0);
1554: }
1556: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt numFields, const PetscInt fields[], PetscSection *subs)
1557: {
1558: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1562: if (!numFields) return(0);
1566: PetscSectionGetNumFields(s, &nF);
1567: if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", numFields, nF);
1568: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1569: PetscSectionSetNumFields(*subs, numFields);
1570: for (f = 0; f < numFields; ++f) {
1571: const char *name = NULL;
1572: PetscInt numComp = 0;
1574: PetscSectionGetFieldName(s, fields[f], &name);
1575: PetscSectionSetFieldName(*subs, f, name);
1576: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1577: PetscSectionSetFieldComponents(*subs, f, numComp);
1578: }
1579: PetscSectionGetChart(s, &pStart, &pEnd);
1580: PetscSectionSetChart(*subs, pStart, pEnd);
1581: for (p = pStart; p < pEnd; ++p) {
1582: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1584: for (f = 0; f < numFields; ++f) {
1585: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1586: PetscSectionSetFieldDof(*subs, p, f, fdof);
1587: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1588: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1589: dof += fdof;
1590: cdof += cfdof;
1591: }
1592: PetscSectionSetDof(*subs, p, dof);
1593: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1594: maxCdof = PetscMax(cdof, maxCdof);
1595: }
1596: PetscSectionSetUp(*subs);
1597: if (maxCdof) {
1598: PetscInt *indices;
1600: PetscMalloc1(maxCdof, &indices);
1601: for (p = pStart; p < pEnd; ++p) {
1602: PetscInt cdof;
1604: PetscSectionGetConstraintDof(*subs, p, &cdof);
1605: if (cdof) {
1606: const PetscInt *oldIndices = NULL;
1607: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1609: for (f = 0; f < numFields; ++f) {
1610: PetscInt oldFoff = 0, oldf;
1612: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1613: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1614: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1615: /* This can be sped up if we assume sorted fields */
1616: for (oldf = 0; oldf < fields[f]; ++oldf) {
1617: PetscInt oldfdof = 0;
1618: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1619: oldFoff += oldfdof;
1620: }
1621: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1622: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1623: numConst += cfdof;
1624: fOff += fdof;
1625: }
1626: PetscSectionSetConstraintIndices(*subs, p, indices);
1627: }
1628: }
1629: PetscFree(indices);
1630: }
1631: return(0);
1632: }
1634: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1635: {
1636: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1640: if (!len) return(0);
1641: for (i = 0; i < len; ++i) {
1642: PetscInt nf, pStarti, pEndi;
1644: PetscSectionGetNumFields(s[i], &nf);
1645: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1646: pStart = PetscMin(pStart, pStarti);
1647: pEnd = PetscMax(pEnd, pEndi);
1648: Nf += nf;
1649: }
1650: PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1651: PetscSectionSetNumFields(*supers, Nf);
1652: for (i = 0, f = 0; i < len; ++i) {
1653: PetscInt nf, fi;
1655: PetscSectionGetNumFields(s[i], &nf);
1656: for (fi = 0; fi < nf; ++fi, ++f) {
1657: const char *name = NULL;
1658: PetscInt numComp = 0;
1660: PetscSectionGetFieldName(s[i], fi, &name);
1661: PetscSectionSetFieldName(*supers, f, name);
1662: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1663: PetscSectionSetFieldComponents(*supers, f, numComp);
1664: }
1665: }
1666: PetscSectionSetChart(*supers, pStart, pEnd);
1667: for (p = pStart; p < pEnd; ++p) {
1668: PetscInt dof = 0, cdof = 0;
1670: for (i = 0, f = 0; i < len; ++i) {
1671: PetscInt nf, fi, pStarti, pEndi;
1672: PetscInt fdof = 0, cfdof = 0;
1674: PetscSectionGetNumFields(s[i], &nf);
1675: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1676: if ((p < pStarti) || (p >= pEndi)) continue;
1677: for (fi = 0; fi < nf; ++fi, ++f) {
1678: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1679: PetscSectionAddFieldDof(*supers, p, f, fdof);
1680: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1681: if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1682: dof += fdof;
1683: cdof += cfdof;
1684: }
1685: }
1686: PetscSectionSetDof(*supers, p, dof);
1687: if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1688: maxCdof = PetscMax(cdof, maxCdof);
1689: }
1690: PetscSectionSetUp(*supers);
1691: if (maxCdof) {
1692: PetscInt *indices;
1694: PetscMalloc1(maxCdof, &indices);
1695: for (p = pStart; p < pEnd; ++p) {
1696: PetscInt cdof;
1698: PetscSectionGetConstraintDof(*supers, p, &cdof);
1699: if (cdof) {
1700: PetscInt dof, numConst = 0, fOff = 0;
1702: for (i = 0, f = 0; i < len; ++i) {
1703: const PetscInt *oldIndices = NULL;
1704: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1706: PetscSectionGetNumFields(s[i], &nf);
1707: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1708: if ((p < pStarti) || (p >= pEndi)) continue;
1709: for (fi = 0; fi < nf; ++fi, ++f) {
1710: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1711: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1712: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1713: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1714: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1715: numConst += cfdof;
1716: }
1717: PetscSectionGetDof(s[i], p, &dof);
1718: fOff += dof;
1719: }
1720: PetscSectionSetConstraintIndices(*supers, p, indices);
1721: }
1722: }
1723: PetscFree(indices);
1724: }
1725: return(0);
1726: }
1728: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1729: {
1730: const PetscInt *points = NULL, *indices = NULL;
1731: PetscInt numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;
1738: PetscSectionGetNumFields(s, &numFields);
1739: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1740: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1741: for (f = 0; f < numFields; ++f) {
1742: const char *name = NULL;
1743: PetscInt numComp = 0;
1745: PetscSectionGetFieldName(s, f, &name);
1746: PetscSectionSetFieldName(*subs, f, name);
1747: PetscSectionGetFieldComponents(s, f, &numComp);
1748: PetscSectionSetFieldComponents(*subs, f, numComp);
1749: }
1750: /* For right now, we do not try to squeeze the subchart */
1751: if (subpointMap) {
1752: ISGetSize(subpointMap, &numSubpoints);
1753: ISGetIndices(subpointMap, &points);
1754: }
1755: PetscSectionGetChart(s, &pStart, &pEnd);
1756: PetscSectionSetChart(*subs, 0, numSubpoints);
1757: for (p = pStart; p < pEnd; ++p) {
1758: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1760: PetscFindInt(p, numSubpoints, points, &subp);
1761: if (subp < 0) continue;
1762: for (f = 0; f < numFields; ++f) {
1763: PetscSectionGetFieldDof(s, p, f, &fdof);
1764: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1765: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1766: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1767: }
1768: PetscSectionGetDof(s, p, &dof);
1769: PetscSectionSetDof(*subs, subp, dof);
1770: PetscSectionGetConstraintDof(s, p, &cdof);
1771: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1772: }
1773: PetscSectionSetUp(*subs);
1774: /* Change offsets to original offsets */
1775: for (p = pStart; p < pEnd; ++p) {
1776: PetscInt off, foff = 0;
1778: PetscFindInt(p, numSubpoints, points, &subp);
1779: if (subp < 0) continue;
1780: for (f = 0; f < numFields; ++f) {
1781: PetscSectionGetFieldOffset(s, p, f, &foff);
1782: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1783: }
1784: PetscSectionGetOffset(s, p, &off);
1785: PetscSectionSetOffset(*subs, subp, off);
1786: }
1787: /* Copy constraint indices */
1788: for (subp = 0; subp < numSubpoints; ++subp) {
1789: PetscInt cdof;
1791: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1792: if (cdof) {
1793: for (f = 0; f < numFields; ++f) {
1794: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1795: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1796: }
1797: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1798: PetscSectionSetConstraintIndices(*subs, subp, indices);
1799: }
1800: }
1801: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1802: return(0);
1803: }
1805: PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1806: {
1807: PetscInt p;
1808: PetscMPIInt rank;
1812: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1813: PetscViewerASCIIPushSynchronized(viewer);
1814: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1815: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1816: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1817: PetscInt b;
1819: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1820: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1821: PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
1822: }
1823: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1824: } else {
1825: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1826: }
1827: }
1828: PetscViewerFlush(viewer);
1829: PetscViewerASCIIPopSynchronized(viewer);
1830: if (s->sym) {
1831: PetscViewerASCIIPushTab(viewer);
1832: PetscSectionSymView(s->sym,viewer);
1833: PetscViewerASCIIPopTab(viewer);
1834: }
1835: return(0);
1836: }
1838: /*@C
1839: PetscSectionView - Views a PetscSection
1841: Collective on PetscSection
1843: Input Parameters:
1844: + s - the PetscSection object to view
1845: - v - the viewer
1847: Level: developer
1849: .seealso PetscSectionCreate(), PetscSectionDestroy()
1850: @*/
1851: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1852: {
1853: PetscBool isascii;
1854: PetscInt f;
1859: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1861: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1862: if (isascii) {
1863: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1864: if (s->numFields) {
1865: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1866: for (f = 0; f < s->numFields; ++f) {
1867: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
1868: PetscSectionView_ASCII(s->field[f], viewer);
1869: }
1870: } else {
1871: PetscSectionView_ASCII(s, viewer);
1872: }
1873: }
1874: return(0);
1875: }
1877: /*@
1878: PetscSectionReset - Frees all section data.
1880: Not collective
1882: Input Parameters:
1883: . s - the PetscSection
1885: Level: developer
1887: .seealso: PetscSection, PetscSectionCreate()
1888: @*/
1889: PetscErrorCode PetscSectionReset(PetscSection s)
1890: {
1891: PetscInt f;
1896: PetscFree(s->numFieldComponents);
1897: for (f = 0; f < s->numFields; ++f) {
1898: PetscSectionDestroy(&s->field[f]);
1899: PetscFree(s->fieldNames[f]);
1900: }
1901: PetscFree(s->fieldNames);
1902: PetscFree(s->field);
1903: PetscSectionDestroy(&s->bc);
1904: PetscFree(s->bcIndices);
1905: PetscFree2(s->atlasDof, s->atlasOff);
1906: PetscSectionDestroy(&s->clSection);
1907: ISDestroy(&s->clPoints);
1908: ISDestroy(&s->perm);
1909: PetscFree(s->clPerm);
1910: PetscFree(s->clInvPerm);
1911: PetscSectionSymDestroy(&s->sym);
1913: s->pStart = -1;
1914: s->pEnd = -1;
1915: s->maxDof = 0;
1916: s->setup = PETSC_FALSE;
1917: s->numFields = 0;
1918: s->clObj = NULL;
1919: return(0);
1920: }
1922: /*@
1923: PetscSectionDestroy - Frees a section object and frees its range if that exists.
1925: Not collective
1927: Input Parameters:
1928: . s - the PetscSection
1930: Level: developer
1932: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
1933: recommended they not be used in user codes unless you really gain something in their use.
1935: .seealso: PetscSection, PetscSectionCreate()
1936: @*/
1937: PetscErrorCode PetscSectionDestroy(PetscSection *s)
1938: {
1942: if (!*s) return(0);
1944: if (--((PetscObject)(*s))->refct > 0) {
1945: *s = NULL;
1946: return(0);
1947: }
1948: PetscSectionReset(*s);
1949: PetscHeaderDestroy(s);
1950: return(0);
1951: }
1953: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
1954: {
1955: const PetscInt p = point - s->pStart;
1959: *values = &baseArray[s->atlasOff[p]];
1960: return(0);
1961: }
1963: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
1964: {
1965: PetscInt *array;
1966: const PetscInt p = point - s->pStart;
1967: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
1968: PetscInt cDim = 0;
1973: PetscSectionGetConstraintDof(s, p, &cDim);
1974: array = &baseArray[s->atlasOff[p]];
1975: if (!cDim) {
1976: if (orientation >= 0) {
1977: const PetscInt dim = s->atlasDof[p];
1978: PetscInt i;
1980: if (mode == INSERT_VALUES) {
1981: for (i = 0; i < dim; ++i) array[i] = values[i];
1982: } else {
1983: for (i = 0; i < dim; ++i) array[i] += values[i];
1984: }
1985: } else {
1986: PetscInt offset = 0;
1987: PetscInt j = -1, field, i;
1989: for (field = 0; field < s->numFields; ++field) {
1990: const PetscInt dim = s->field[field]->atlasDof[p];
1992: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
1993: offset += dim;
1994: }
1995: }
1996: } else {
1997: if (orientation >= 0) {
1998: const PetscInt dim = s->atlasDof[p];
1999: PetscInt cInd = 0, i;
2000: const PetscInt *cDof;
2002: PetscSectionGetConstraintIndices(s, point, &cDof);
2003: if (mode == INSERT_VALUES) {
2004: for (i = 0; i < dim; ++i) {
2005: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2006: array[i] = values[i];
2007: }
2008: } else {
2009: for (i = 0; i < dim; ++i) {
2010: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2011: array[i] += values[i];
2012: }
2013: }
2014: } else {
2015: const PetscInt *cDof;
2016: PetscInt offset = 0;
2017: PetscInt cOffset = 0;
2018: PetscInt j = 0, field;
2020: PetscSectionGetConstraintIndices(s, point, &cDof);
2021: for (field = 0; field < s->numFields; ++field) {
2022: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2023: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2024: const PetscInt sDim = dim - tDim;
2025: PetscInt cInd = 0, i,k;
2027: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2028: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2029: array[j] = values[k];
2030: }
2031: offset += dim;
2032: cOffset += dim - tDim;
2033: }
2034: }
2035: }
2036: return(0);
2037: }
2039: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2040: {
2044: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2045: return(0);
2046: }
2048: /*@C
2049: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2051: Input Parameters:
2052: + s - The PetscSection
2053: - point - The point
2055: Output Parameter:
2056: . indices - The constrained dofs
2058: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
2060: Level: advanced
2062: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2063: @*/
2064: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2065: {
2070: if (s->bc) {
2071: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2072: } else *indices = NULL;
2073: return(0);
2074: }
2076: /*@C
2077: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2079: Input Parameters:
2080: + s - The PetscSection
2081: . point - The point
2082: - indices - The constrained dofs
2084: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
2086: Level: advanced
2088: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2089: @*/
2090: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2091: {
2096: if (s->bc) {
2097: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2098: }
2099: return(0);
2100: }
2102: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2103: {
2108: 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);
2109: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2110: return(0);
2111: }
2113: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2114: {
2119: 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);
2120: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2121: return(0);
2122: }
2124: /*@
2125: PetscSectionPermute - Reorder the section according to the input point permutation
2127: Collective on PetscSection
2129: Input Parameter:
2130: + section - The PetscSection object
2131: - perm - The point permutation, old point p becomes new point perm[p]
2133: Output Parameter:
2134: . sectionNew - The permuted PetscSection
2136: Level: intermediate
2138: .keywords: mesh
2139: .seealso: MatPermute()
2140: @*/
2141: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2142: {
2143: PetscSection s = section, sNew;
2144: const PetscInt *perm;
2145: PetscInt numFields, f, numPoints, pStart, pEnd, p;
2146: PetscErrorCode ierr;
2152: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2153: PetscSectionGetNumFields(s, &numFields);
2154: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2155: for (f = 0; f < numFields; ++f) {
2156: const char *name;
2157: PetscInt numComp;
2159: PetscSectionGetFieldName(s, f, &name);
2160: PetscSectionSetFieldName(sNew, f, name);
2161: PetscSectionGetFieldComponents(s, f, &numComp);
2162: PetscSectionSetFieldComponents(sNew, f, numComp);
2163: }
2164: ISGetLocalSize(permutation, &numPoints);
2165: ISGetIndices(permutation, &perm);
2166: PetscSectionGetChart(s, &pStart, &pEnd);
2167: PetscSectionSetChart(sNew, pStart, pEnd);
2168: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2169: for (p = pStart; p < pEnd; ++p) {
2170: PetscInt dof, cdof;
2172: PetscSectionGetDof(s, p, &dof);
2173: PetscSectionSetDof(sNew, perm[p], dof);
2174: PetscSectionGetConstraintDof(s, p, &cdof);
2175: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2176: for (f = 0; f < numFields; ++f) {
2177: PetscSectionGetFieldDof(s, p, f, &dof);
2178: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2179: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2180: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2181: }
2182: }
2183: PetscSectionSetUp(sNew);
2184: for (p = pStart; p < pEnd; ++p) {
2185: const PetscInt *cind;
2186: PetscInt cdof;
2188: PetscSectionGetConstraintDof(s, p, &cdof);
2189: if (cdof) {
2190: PetscSectionGetConstraintIndices(s, p, &cind);
2191: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2192: }
2193: for (f = 0; f < numFields; ++f) {
2194: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2195: if (cdof) {
2196: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2197: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2198: }
2199: }
2200: }
2201: ISRestoreIndices(permutation, &perm);
2202: *sectionNew = sNew;
2203: return(0);
2204: }
2206: /*@C
2207: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
2209: Collective
2211: Input Parameters:
2212: + sf - The SF
2213: - rootSection - Section defined on root space
2215: Output Parameters:
2216: + remoteOffsets - root offsets in leaf storage, or NULL
2217: - leafSection - Section defined on the leaf space
2219: Level: intermediate
2221: .seealso: PetscSFCreate()
2222: @*/
2223: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2224: {
2225: PetscSF embedSF;
2226: const PetscInt *ilocal, *indices;
2227: IS selected;
2228: PetscInt numFields, nroots, nleaves, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, i, f;
2232: PetscSectionGetNumFields(rootSection, &numFields);
2233: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2234: for (f = 0; f < numFields; ++f) {
2235: PetscInt numComp = 0;
2236: PetscSectionGetFieldComponents(rootSection, f, &numComp);
2237: PetscSectionSetFieldComponents(leafSection, f, numComp);
2238: }
2239: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2240: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2241: rpEnd = PetscMin(rpEnd,nroots);
2242: rpEnd = PetscMax(rpStart,rpEnd);
2243: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2244: ISGetIndices(selected, &indices);
2245: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2246: ISRestoreIndices(selected, &indices);
2247: ISDestroy(&selected);
2248: PetscSFGetGraph(embedSF, NULL, &nleaves, &ilocal, NULL);
2249: if (nleaves && ilocal) {
2250: for (i = 0; i < nleaves; ++i) {
2251: lpStart = PetscMin(lpStart, ilocal[i]);
2252: lpEnd = PetscMax(lpEnd, ilocal[i]);
2253: }
2254: ++lpEnd;
2255: } else {
2256: lpStart = 0;
2257: lpEnd = nleaves;
2258: }
2259: PetscSectionSetChart(leafSection, lpStart, lpEnd);
2260: /* Could fuse these at the cost of a copy and extra allocation */
2261: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2262: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2263: for (f = 0; f < numFields; ++f) {
2264: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2265: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2266: }
2267: if (remoteOffsets) {
2268: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2269: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2270: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2271: }
2272: PetscSFDestroy(&embedSF);
2273: PetscSectionSetUp(leafSection);
2274: return(0);
2275: }
2277: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2278: {
2279: PetscSF embedSF;
2280: const PetscInt *indices;
2281: IS selected;
2282: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2283: PetscErrorCode ierr;
2286: *remoteOffsets = NULL;
2287: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2288: if (numRoots < 0) return(0);
2289: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2290: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2291: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2292: ISGetIndices(selected, &indices);
2293: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2294: ISRestoreIndices(selected, &indices);
2295: ISDestroy(&selected);
2296: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2297: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2298: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2299: PetscSFDestroy(&embedSF);
2300: return(0);
2301: }
2303: /*@C
2304: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2306: Input Parameters:
2307: + sf - The SF
2308: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2309: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2310: - leafSection - Data layout of local points for incoming data (this is the distributed section)
2312: Output Parameters:
2313: - sectionSF - The new SF
2315: Note: Either rootSection or remoteOffsets can be specified
2317: Level: intermediate
2319: .seealso: PetscSFCreate()
2320: @*/
2321: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2322: {
2323: MPI_Comm comm;
2324: const PetscInt *localPoints;
2325: const PetscSFNode *remotePoints;
2326: PetscInt lpStart, lpEnd;
2327: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2328: PetscInt *localIndices;
2329: PetscSFNode *remoteIndices;
2330: PetscInt i, ind;
2331: PetscErrorCode ierr;
2339: PetscObjectGetComm((PetscObject)sf,&comm);
2340: PetscSFCreate(comm, sectionSF);
2341: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2342: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2343: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2344: if (numRoots < 0) return(0);
2345: for (i = 0; i < numPoints; ++i) {
2346: PetscInt localPoint = localPoints ? localPoints[i] : i;
2347: PetscInt dof;
2349: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2350: PetscSectionGetDof(leafSection, localPoint, &dof);
2351: numIndices += dof;
2352: }
2353: }
2354: PetscMalloc1(numIndices, &localIndices);
2355: PetscMalloc1(numIndices, &remoteIndices);
2356: /* Create new index graph */
2357: for (i = 0, ind = 0; i < numPoints; ++i) {
2358: PetscInt localPoint = localPoints ? localPoints[i] : i;
2359: PetscInt rank = remotePoints[i].rank;
2361: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2362: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2363: PetscInt loff, dof, d;
2365: PetscSectionGetOffset(leafSection, localPoint, &loff);
2366: PetscSectionGetDof(leafSection, localPoint, &dof);
2367: for (d = 0; d < dof; ++d, ++ind) {
2368: localIndices[ind] = loff+d;
2369: remoteIndices[ind].rank = rank;
2370: remoteIndices[ind].index = remoteOffset+d;
2371: }
2372: }
2373: }
2374: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2375: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2376: return(0);
2377: }
2379: /*@
2380: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2382: Input Parameters:
2383: + section - The PetscSection
2384: . obj - A PetscObject which serves as the key for this index
2385: . clSection - Section giving the size of the closure of each point
2386: - clPoints - IS giving the points in each closure
2388: Note: We compress out closure points with no dofs in this section
2390: Level: intermediate
2392: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2393: @*/
2394: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2395: {
2399: if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2400: section->clObj = obj;
2401: PetscSectionDestroy(§ion->clSection);
2402: ISDestroy(§ion->clPoints);
2403: section->clSection = clSection;
2404: section->clPoints = clPoints;
2405: return(0);
2406: }
2408: /*@
2409: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2411: Input Parameters:
2412: + section - The PetscSection
2413: - obj - A PetscObject which serves as the key for this index
2415: Output Parameters:
2416: + clSection - Section giving the size of the closure of each point
2417: - clPoints - IS giving the points in each closure
2419: Note: We compress out closure points with no dofs in this section
2421: Level: intermediate
2423: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2424: @*/
2425: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2426: {
2428: if (section->clObj == obj) {
2429: if (clSection) *clSection = section->clSection;
2430: if (clPoints) *clPoints = section->clPoints;
2431: } else {
2432: if (clSection) *clSection = NULL;
2433: if (clPoints) *clPoints = NULL;
2434: }
2435: return(0);
2436: }
2438: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2439: {
2440: PetscInt i;
2444: if (section->clObj != obj) {
2445: PetscSectionDestroy(§ion->clSection);
2446: ISDestroy(§ion->clPoints);
2447: }
2448: section->clObj = obj;
2449: PetscFree(section->clPerm);
2450: PetscFree(section->clInvPerm);
2451: section->clSize = clSize;
2452: if (mode == PETSC_COPY_VALUES) {
2453: PetscMalloc1(clSize, §ion->clPerm);
2454: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2455: PetscMemcpy(section->clPerm, clPerm, clSize*sizeof(PetscInt));
2456: } else if (mode == PETSC_OWN_POINTER) {
2457: section->clPerm = clPerm;
2458: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2459: PetscMalloc1(clSize, §ion->clInvPerm);
2460: for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2461: return(0);
2462: }
2464: /*@
2465: PetscSectionSetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2467: Input Parameters:
2468: + section - The PetscSection
2469: . obj - A PetscObject which serves as the key for this index
2470: - perm - Permutation of the cell dof closure
2472: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2473: other points (like faces).
2475: Level: intermediate
2477: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2478: @*/
2479: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2480: {
2481: const PetscInt *clPerm = NULL;
2482: PetscInt clSize = 0;
2483: PetscErrorCode ierr;
2486: if (perm) {
2487: ISGetLocalSize(perm, &clSize);
2488: ISGetIndices(perm, &clPerm);
2489: }
2490: PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2491: if (perm) {ISRestoreIndices(perm, &clPerm);}
2492: return(0);
2493: }
2495: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2496: {
2498: if (section->clObj == obj) {
2499: if (size) *size = section->clSize;
2500: if (perm) *perm = section->clPerm;
2501: } else {
2502: if (size) *size = 0;
2503: if (perm) *perm = NULL;
2504: }
2505: return(0);
2506: }
2508: /*@
2509: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2511: Input Parameters:
2512: + section - The PetscSection
2513: - obj - A PetscObject which serves as the key for this index
2515: Output Parameter:
2516: . perm - The dof closure permutation
2518: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2519: other points (like faces).
2521: The user must destroy the IS that is returned.
2523: Level: intermediate
2525: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2526: @*/
2527: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2528: {
2529: const PetscInt *clPerm;
2530: PetscInt clSize;
2531: PetscErrorCode ierr;
2534: PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2535: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2536: return(0);
2537: }
2539: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2540: {
2542: if (section->clObj == obj) {
2543: if (size) *size = section->clSize;
2544: if (perm) *perm = section->clInvPerm;
2545: } else {
2546: if (size) *size = 0;
2547: if (perm) *perm = NULL;
2548: }
2549: return(0);
2550: }
2552: /*@
2553: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2555: Input Parameters:
2556: + section - The PetscSection
2557: - obj - A PetscObject which serves as the key for this index
2559: Output Parameters:
2560: + size - The dof closure size
2561: - perm - The dof closure permutation
2563: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2564: other points (like faces).
2566: The user must destroy the IS that is returned.
2568: Level: intermediate
2570: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2571: @*/
2572: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2573: {
2574: const PetscInt *clPerm;
2575: PetscInt clSize;
2576: PetscErrorCode ierr;
2579: PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2580: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2581: return(0);
2582: }
2584: /*@
2585: PetscSectionGetField - Get the subsection associated with a single field
2587: Input Parameters:
2588: + s - The PetscSection
2589: - field - The field number
2591: Output Parameter:
2592: . subs - The subsection for the given field
2594: Level: intermediate
2596: .seealso: PetscSectionSetNumFields()
2597: @*/
2598: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2599: {
2603: 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);
2604: *subs = s->field[field];
2605: return(0);
2606: }
2608: PetscClassId PETSC_SECTION_SYM_CLASSID;
2609: PetscFunctionList PetscSectionSymList = NULL;
2611: /*@
2612: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
2614: Collective on MPI_Comm
2616: Input Parameter:
2617: . comm - the MPI communicator
2619: Output Parameter:
2620: . sym - pointer to the new set of symmetries
2622: Level: developer
2624: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2625: @*/
2626: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2627: {
2632: ISInitializePackage();
2633: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2634: return(0);
2635: }
2637: /*@C
2638: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
2640: Collective on PetscSectionSym
2642: Input Parameters:
2643: + sym - The section symmetry object
2644: - method - The name of the section symmetry type
2646: Level: developer
2648: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2649: @*/
2650: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2651: {
2652: PetscErrorCode (*r)(PetscSectionSym);
2653: PetscBool match;
2658: PetscObjectTypeCompare((PetscObject) sym, method, &match);
2659: if (match) return(0);
2661: PetscFunctionListFind(PetscSectionSymList,method,&r);
2662: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2663: if (sym->ops->destroy) {
2664: (*sym->ops->destroy)(sym);
2665: sym->ops->destroy = NULL;
2666: }
2667: (*r)(sym);
2668: PetscObjectChangeTypeName((PetscObject)sym,method);
2669: return(0);
2670: }
2673: /*@C
2674: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
2676: Not Collective
2678: Input Parameter:
2679: . sym - The section symmetry
2681: Output Parameter:
2682: . type - The index set type name
2684: Level: developer
2686: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2687: @*/
2688: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2689: {
2693: *type = ((PetscObject)sym)->type_name;
2694: return(0);
2695: }
2697: /*@C
2698: PetscSectionSymRegister - Adds a new section symmetry implementation
2700: Not Collective
2702: Input Parameters:
2703: + name - The name of a new user-defined creation routine
2704: - create_func - The creation routine itself
2706: Notes:
2707: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
2709: Level: developer
2711: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2712: @*/
2713: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2714: {
2718: ISInitializePackage();
2719: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2720: return(0);
2721: }
2723: /*@
2724: PetscSectionSymDestroy - Destroys a section symmetry.
2726: Collective on PetscSectionSym
2728: Input Parameters:
2729: . sym - the section symmetry
2731: Level: developer
2733: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2734: @*/
2735: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2736: {
2737: SymWorkLink link,next;
2741: if (!*sym) return(0);
2743: if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
2744: if ((*sym)->ops->destroy) {
2745: (*(*sym)->ops->destroy)(*sym);
2746: }
2747: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
2748: for (link=(*sym)->workin; link; link=next) {
2749: next = link->next;
2750: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2751: PetscFree(link);
2752: }
2753: (*sym)->workin = NULL;
2754: PetscHeaderDestroy(sym);
2755: return(0);
2756: }
2758: /*@C
2759: PetscSectionSymView - Displays a section symmetry
2761: Collective on PetscSectionSym
2763: Input Parameters:
2764: + sym - the index set
2765: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
2767: Level: developer
2769: .seealso: PetscViewerASCIIOpen()
2770: @*/
2771: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
2772: {
2777: if (!viewer) {
2778: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
2779: }
2782: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
2783: if (sym->ops->view) {
2784: (*sym->ops->view)(sym,viewer);
2785: }
2786: return(0);
2787: }
2789: /*@
2790: PetscSectionSetSym - Set the symmetries for the data referred to by the section
2792: Collective on PetscSection
2794: Input Parameters:
2795: + section - the section describing data layout
2796: - sym - the symmetry describing the affect of orientation on the access of the data
2798: Level: developer
2800: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
2801: @*/
2802: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
2803: {
2808: PetscSectionSymDestroy(&(section->sym));
2809: if (sym) {
2812: PetscObjectReference((PetscObject) sym);
2813: }
2814: section->sym = sym;
2815: return(0);
2816: }
2818: /*@
2819: PetscSectionGetSym - Get the symmetries for the data referred to by the section
2821: Not collective
2823: Input Parameters:
2824: . section - the section describing data layout
2826: Output Parameters:
2827: . sym - the symmetry describing the affect of orientation on the access of the data
2829: Level: developer
2831: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
2832: @*/
2833: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
2834: {
2837: *sym = section->sym;
2838: return(0);
2839: }
2841: /*@
2842: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
2844: Collective on PetscSection
2846: Input Parameters:
2847: + section - the section describing data layout
2848: . field - the field number
2849: - sym - the symmetry describing the affect of orientation on the access of the data
2851: Level: developer
2853: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
2854: @*/
2855: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
2856: {
2861: 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);
2862: PetscSectionSetSym(section->field[field],sym);
2863: return(0);
2864: }
2866: /*@
2867: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
2869: Collective on PetscSection
2871: Input Parameters:
2872: + section - the section describing data layout
2873: - field - the field number
2875: Output Parameters:
2876: . sym - the symmetry describing the affect of orientation on the access of the data
2878: Level: developer
2880: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
2881: @*/
2882: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
2883: {
2886: 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);
2887: *sym = section->field[field]->sym;
2888: return(0);
2889: }
2891: /*@C
2892: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
2894: Not collective
2896: Input Parameters:
2897: + section - the section
2898: . numPoints - the number of points
2899: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
2900: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
2901: context, see DMPlexGetConeOrientation()).
2903: Output Parameter:
2904: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
2905: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
2906: identity).
2908: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
2909: .vb
2910: const PetscInt **perms;
2911: const PetscScalar **rots;
2912: PetscInt lOffset;
2914: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2915: for (i = 0, lOffset = 0; i < numPoints; i++) {
2916: PetscInt point = points[2*i], dof, sOffset;
2917: const PetscInt *perm = perms ? perms[i] : NULL;
2918: const PetscScalar *rot = rots ? rots[i] : NULL;
2920: PetscSectionGetDof(section,point,&dof);
2921: PetscSectionGetOffset(section,point,&sOffset);
2923: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
2924: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
2925: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
2926: lOffset += dof;
2927: }
2928: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2929: .ve
2931: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
2932: .vb
2933: const PetscInt **perms;
2934: const PetscScalar **rots;
2935: PetscInt lOffset;
2937: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
2938: for (i = 0, lOffset = 0; i < numPoints; i++) {
2939: PetscInt point = points[2*i], dof, sOffset;
2940: const PetscInt *perm = perms ? perms[i] : NULL;
2941: const PetscScalar *rot = rots ? rots[i] : NULL;
2943: PetscSectionGetDof(section,point,&dof);
2944: PetscSectionGetOffset(section,point,&sOff);
2946: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
2947: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
2948: offset += dof;
2949: }
2950: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
2951: .ve
2953: Level: developer
2955: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
2956: @*/
2957: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
2958: {
2959: PetscSectionSym sym;
2960: PetscErrorCode ierr;
2965: if (perms) *perms = NULL;
2966: if (rots) *rots = NULL;
2967: sym = section->sym;
2968: if (sym && (perms || rots)) {
2969: SymWorkLink link;
2971: if (sym->workin) {
2972: link = sym->workin;
2973: sym->workin = sym->workin->next;
2974: } else {
2975: PetscNewLog(sym,&link);
2976: }
2977: if (numPoints > link->numPoints) {
2978: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
2979: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
2980: link->numPoints = numPoints;
2981: }
2982: link->next = sym->workout;
2983: sym->workout = link;
2984: PetscMemzero((void *) link->perms,numPoints * sizeof(const PetscInt *));
2985: PetscMemzero((void *) link->rots,numPoints * sizeof(const PetscScalar *));
2986: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
2987: if (perms) *perms = link->perms;
2988: if (rots) *rots = link->rots;
2989: }
2990: return(0);
2991: }
2993: /*@C
2994: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
2996: Not collective
2998: Input Parameters:
2999: + section - the section
3000: . numPoints - the number of points
3001: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3002: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3003: context, see DMPlexGetConeOrientation()).
3005: Output Parameter:
3006: + perms - The permutations for the given orientations: set to NULL at conclusion
3007: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3009: Level: developer
3011: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3012: @*/
3013: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3014: {
3015: PetscSectionSym sym;
3019: sym = section->sym;
3020: if (sym && (perms || rots)) {
3021: SymWorkLink *p,link;
3023: for (p=&sym->workout; (link=*p); p=&link->next) {
3024: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3025: *p = link->next;
3026: link->next = sym->workin;
3027: sym->workin = link;
3028: if (perms) *perms = NULL;
3029: if (rots) *rots = NULL;
3030: return(0);
3031: }
3032: }
3033: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3034: }
3035: return(0);
3036: }
3038: /*@C
3039: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
3041: Not collective
3043: Input Parameters:
3044: + section - the section
3045: . field - the field of the section
3046: . numPoints - the number of points
3047: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3048: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3049: context, see DMPlexGetConeOrientation()).
3051: Output Parameter:
3052: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3053: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3054: identity).
3056: Level: developer
3058: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3059: @*/
3060: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3061: {
3066: 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);
3067: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3068: return(0);
3069: }
3071: /*@C
3072: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
3074: Not collective
3076: Input Parameters:
3077: + section - the section
3078: . field - the field number
3079: . numPoints - the number of points
3080: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3081: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3082: context, see DMPlexGetConeOrientation()).
3084: Output Parameter:
3085: + perms - The permutations for the given orientations: set to NULL at conclusion
3086: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3088: Level: developer
3090: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3091: @*/
3092: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3093: {
3098: 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);
3099: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3100: return(0);
3101: }
3103: /*@
3104: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3106: Not collective
3108: Input Parameter:
3109: . s - the global PetscSection
3111: Output Parameters:
3112: . flg - the flag
3114: Level: developer
3116: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3117: @*/
3118: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3119: {
3122: *flg = s->useFieldOff;
3123: return(0);
3124: }
3126: /*@
3127: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3129: Not collective
3131: Input Parameters:
3132: + s - the global PetscSection
3133: - flg - the flag
3135: Level: developer
3137: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3138: @*/
3139: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3140: {
3143: s->useFieldOff = flg;
3144: return(0);
3145: }
3147: #define PetscSectionExpandPoints_Loop(TYPE) \
3148: { \
3149: PetscInt i, n, o0, o1, size; \
3150: TYPE *a0 = (TYPE*)origArray, *a1; \
3151: PetscSectionGetStorageSize(s, &size); \
3152: PetscMalloc1(size, &a1); \
3153: for (i=0; i<npoints; i++) { \
3154: PetscSectionGetOffset(origSection, points_[i], &o0); \
3155: PetscSectionGetOffset(s, i, &o1); \
3156: PetscSectionGetDof(s, i, &n); \
3157: PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3158: } \
3159: *newArray = (void*)a1; \
3160: }
3162: /*@
3163: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3165: Not collective
3167: Input Parameters:
3168: + origSection - the PetscSection describing the layout of the array
3169: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3170: . origArray - the array; its size must be equal to the storage size of origSection
3171: - points - IS with points to extract; its indices must lie in the chart of origSection
3173: Output Parameters:
3174: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3175: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3177: Level: developer
3179: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3180: @*/
3181: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3182: {
3183: PetscSection s;
3184: const PetscInt *points_;
3185: PetscInt i, n, npoints, pStart, pEnd;
3186: PetscMPIInt unitsize;
3187: PetscErrorCode ierr;
3195: MPI_Type_size(dataType, &unitsize);
3196: ISGetLocalSize(points, &npoints);
3197: ISGetIndices(points, &points_);
3198: PetscSectionGetChart(origSection, &pStart, &pEnd);
3199: PetscSectionCreate(PETSC_COMM_SELF, &s);
3200: PetscSectionSetChart(s, 0, npoints);
3201: for (i=0; i<npoints; i++) {
3202: 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);
3203: PetscSectionGetDof(origSection, points_[i], &n);
3204: PetscSectionSetDof(s, i, n);
3205: }
3206: PetscSectionSetUp(s);
3207: if (newArray) {
3208: if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);}
3209: else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);}
3210: else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);}
3211: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3212: }
3213: if (newSection) {
3214: *newSection = s;
3215: } else {
3216: PetscSectionDestroy(&s);
3217: }
3218: return(0);
3219: }