Actual source code: section.c
petsc-3.13.6 2020-09-29
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsection.h>
7: #include <petscsf.h>
8: #include <petscviewer.h>
10: PetscClassId PETSC_SECTION_CLASSID;
12: /*@
13: PetscSectionCreate - Allocates PetscSection space and sets the map contents to the default.
15: Collective
17: Input Parameters:
18: + comm - the MPI communicator
19: - s - pointer to the section
21: Level: beginner
23: Notes:
24: Typical calling sequence
25: $ PetscSectionCreate(MPI_Comm,PetscSection *);
26: $ PetscSectionSetNumFields(PetscSection, numFields);
27: $ PetscSectionSetChart(PetscSection,low,high);
28: $ PetscSectionSetDof(PetscSection,point,numdof);
29: $ PetscSectionSetUp(PetscSection);
30: $ PetscSectionGetOffset(PetscSection,point,PetscInt *);
31: $ PetscSectionDestroy(PetscSection);
33: The PetscSection object and methods are intended to be used in the PETSc Vec and Mat implementions; it is
34: recommended they not be used in user codes unless you really gain something in their use.
36: .seealso: PetscSection, PetscSectionDestroy()
37: @*/
38: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
39: {
44: ISInitializePackage();
46: PetscHeaderCreate(*s,PETSC_SECTION_CLASSID,"PetscSection","Section","IS",comm,PetscSectionDestroy,PetscSectionView);
48: (*s)->pStart = -1;
49: (*s)->pEnd = -1;
50: (*s)->perm = NULL;
51: (*s)->pointMajor = PETSC_TRUE;
52: (*s)->maxDof = 0;
53: (*s)->atlasDof = NULL;
54: (*s)->atlasOff = NULL;
55: (*s)->bc = NULL;
56: (*s)->bcIndices = NULL;
57: (*s)->setup = PETSC_FALSE;
58: (*s)->numFields = 0;
59: (*s)->fieldNames = NULL;
60: (*s)->field = NULL;
61: (*s)->useFieldOff = PETSC_FALSE;
62: (*s)->compNames = NULL;
63: (*s)->clObj = NULL;
64: (*s)->clSection = NULL;
65: (*s)->clPoints = NULL;
66: (*s)->clSize = 0;
67: (*s)->clPerm = NULL;
68: (*s)->clInvPerm = NULL;
69: return(0);
70: }
72: /*@
73: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
75: Collective
77: Input Parameter:
78: . section - the PetscSection
80: Output Parameter:
81: . newSection - the copy
83: Level: intermediate
85: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
86: @*/
87: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
88: {
89: PetscSectionSym sym;
90: IS perm;
91: PetscInt numFields, f, c, pStart, pEnd, p;
92: PetscErrorCode ierr;
97: PetscSectionReset(newSection);
98: PetscSectionGetNumFields(section, &numFields);
99: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
100: for (f = 0; f < numFields; ++f) {
101: const char *fieldName = NULL, *compName = NULL;
102: PetscInt numComp = 0;
104: PetscSectionGetFieldName(section, f, &fieldName);
105: PetscSectionSetFieldName(newSection, f, fieldName);
106: PetscSectionGetFieldComponents(section, f, &numComp);
107: PetscSectionSetFieldComponents(newSection, f, numComp);
108: for (c = 0; c < numComp; ++c) {
109: PetscSectionGetComponentName(section, f, c, &compName);
110: PetscSectionSetComponentName(newSection, f, c, compName);
111: }
112: PetscSectionGetFieldSym(section, f, &sym);
113: PetscSectionSetFieldSym(newSection, f, sym);
114: }
115: PetscSectionGetChart(section, &pStart, &pEnd);
116: PetscSectionSetChart(newSection, pStart, pEnd);
117: PetscSectionGetPermutation(section, &perm);
118: PetscSectionSetPermutation(newSection, perm);
119: PetscSectionGetSym(section, &sym);
120: PetscSectionSetSym(newSection, sym);
121: for (p = pStart; p < pEnd; ++p) {
122: PetscInt dof, cdof, fcdof = 0;
124: PetscSectionGetDof(section, p, &dof);
125: PetscSectionSetDof(newSection, p, dof);
126: PetscSectionGetConstraintDof(section, p, &cdof);
127: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
128: for (f = 0; f < numFields; ++f) {
129: PetscSectionGetFieldDof(section, p, f, &dof);
130: PetscSectionSetFieldDof(newSection, p, f, dof);
131: if (cdof) {
132: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
133: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
134: }
135: }
136: }
137: PetscSectionSetUp(newSection);
138: for (p = pStart; p < pEnd; ++p) {
139: PetscInt off, cdof, fcdof = 0;
140: const PetscInt *cInd;
142: /* Must set offsets in case they do not agree with the prefix sums */
143: PetscSectionGetOffset(section, p, &off);
144: PetscSectionSetOffset(newSection, p, off);
145: PetscSectionGetConstraintDof(section, p, &cdof);
146: if (cdof) {
147: PetscSectionGetConstraintIndices(section, p, &cInd);
148: PetscSectionSetConstraintIndices(newSection, p, cInd);
149: for (f = 0; f < numFields; ++f) {
150: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
151: if (fcdof) {
152: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
153: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
154: }
155: }
156: }
157: }
158: return(0);
159: }
161: /*@
162: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
164: Collective
166: Input Parameter:
167: . section - the PetscSection
169: Output Parameter:
170: . newSection - the copy
172: Level: beginner
174: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
175: @*/
176: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
177: {
183: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
184: PetscSectionCopy(section, *newSection);
185: return(0);
186: }
188: /*@
189: PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database
191: Collective on PetscSection
193: Input Parameter:
194: . section - the PetscSection
196: Options Database:
197: . -petscsection_point_major the dof order
199: Level: intermediate
201: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
202: @*/
203: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
204: {
209: PetscObjectOptionsBegin((PetscObject) s);
210: PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
211: /* process any options handlers added with PetscObjectAddOptionsHandler() */
212: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
213: PetscOptionsEnd();
214: PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
215: return(0);
216: }
218: /*@
219: PetscSectionCompare - Compares two sections
221: Collective on PetscSection
223: Input Parameters:
224: + s1 - the first PetscSection
225: - s2 - the second PetscSection
227: Output Parameter:
228: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
230: Level: intermediate
232: Notes:
233: Field names are disregarded.
235: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
236: @*/
237: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
238: {
239: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
240: const PetscInt *idx1, *idx2;
241: IS perm1, perm2;
242: PetscBool flg;
243: PetscMPIInt mflg;
250: flg = PETSC_FALSE;
252: MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
253: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
254: *congruent = PETSC_FALSE;
255: return(0);
256: }
258: PetscSectionGetChart(s1, &pStart, &pEnd);
259: PetscSectionGetChart(s2, &n1, &n2);
260: if (pStart != n1 || pEnd != n2) goto not_congruent;
262: PetscSectionGetPermutation(s1, &perm1);
263: PetscSectionGetPermutation(s2, &perm2);
264: if (perm1 && perm2) {
265: ISEqual(perm1, perm2, congruent);
266: if (!(*congruent)) goto not_congruent;
267: } else if (perm1 != perm2) goto not_congruent;
269: for (p = pStart; p < pEnd; ++p) {
270: PetscSectionGetOffset(s1, p, &n1);
271: PetscSectionGetOffset(s2, p, &n2);
272: if (n1 != n2) goto not_congruent;
274: PetscSectionGetDof(s1, p, &n1);
275: PetscSectionGetDof(s2, p, &n2);
276: if (n1 != n2) goto not_congruent;
278: PetscSectionGetConstraintDof(s1, p, &ncdof);
279: PetscSectionGetConstraintDof(s2, p, &n2);
280: if (ncdof != n2) goto not_congruent;
282: PetscSectionGetConstraintIndices(s1, p, &idx1);
283: PetscSectionGetConstraintIndices(s2, p, &idx2);
284: PetscArraycmp(idx1, idx2, ncdof, congruent);
285: if (!(*congruent)) goto not_congruent;
286: }
288: PetscSectionGetNumFields(s1, &nfields);
289: PetscSectionGetNumFields(s2, &n2);
290: if (nfields != n2) goto not_congruent;
292: for (f = 0; f < nfields; ++f) {
293: PetscSectionGetFieldComponents(s1, f, &n1);
294: PetscSectionGetFieldComponents(s2, f, &n2);
295: if (n1 != n2) goto not_congruent;
297: for (p = pStart; p < pEnd; ++p) {
298: PetscSectionGetFieldOffset(s1, p, f, &n1);
299: PetscSectionGetFieldOffset(s2, p, f, &n2);
300: if (n1 != n2) goto not_congruent;
302: PetscSectionGetFieldDof(s1, p, f, &n1);
303: PetscSectionGetFieldDof(s2, p, f, &n2);
304: if (n1 != n2) goto not_congruent;
306: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
307: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
308: if (nfcdof != n2) goto not_congruent;
310: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
311: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
312: PetscArraycmp(idx1, idx2, nfcdof, congruent);
313: if (!(*congruent)) goto not_congruent;
314: }
315: }
317: flg = PETSC_TRUE;
318: not_congruent:
319: MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
320: return(0);
321: }
323: /*@
324: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
326: Not collective
328: Input Parameter:
329: . s - the PetscSection
331: Output Parameter:
332: . numFields - the number of fields defined, or 0 if none were defined
334: Level: intermediate
336: .seealso: PetscSectionSetNumFields()
337: @*/
338: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
339: {
343: *numFields = s->numFields;
344: return(0);
345: }
347: /*@
348: PetscSectionSetNumFields - Sets the number of fields.
350: Not collective
352: Input Parameters:
353: + s - the PetscSection
354: - numFields - the number of fields
356: Level: intermediate
358: .seealso: PetscSectionGetNumFields()
359: @*/
360: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
361: {
362: PetscInt f;
367: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %D must be positive", numFields);
368: PetscSectionReset(s);
370: s->numFields = numFields;
371: PetscMalloc1(s->numFields, &s->numFieldComponents);
372: PetscMalloc1(s->numFields, &s->fieldNames);
373: PetscMalloc1(s->numFields, &s->compNames);
374: PetscMalloc1(s->numFields, &s->field);
375: for (f = 0; f < s->numFields; ++f) {
376: char name[64];
378: s->numFieldComponents[f] = 1;
380: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
381: PetscSNPrintf(name, 64, "Field_%D", f);
382: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
383: PetscSNPrintf(name, 64, "Component_0");
384: PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]);
385: PetscStrallocpy(name, (char **) &s->compNames[f][0]);
386: }
387: return(0);
388: }
390: /*@C
391: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
393: Not Collective
395: Input Parameters:
396: + s - the PetscSection
397: - field - the field number
399: Output Parameter:
400: . fieldName - the field name
402: Level: intermediate
404: .seealso: PetscSectionSetFieldName()
405: @*/
406: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
407: {
411: 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);
412: *fieldName = s->fieldNames[field];
413: return(0);
414: }
416: /*@C
417: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
419: Not Collective
421: Input Parameters:
422: + s - the PetscSection
423: . field - the field number
424: - fieldName - the field name
426: Level: intermediate
428: .seealso: PetscSectionGetFieldName()
429: @*/
430: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
431: {
437: 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);
438: PetscFree(s->fieldNames[field]);
439: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
440: return(0);
441: }
443: /*@C
444: PetscSectionGetComponentName - Gets the name of a field component in the PetscSection
446: Not Collective
448: Input Parameters:
449: + s - the PetscSection
450: . field - the field number
451: . comp - the component number
452: - compName - the component name
454: Level: intermediate
456: .seealso: PetscSectionSetComponentName()
457: @*/
458: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
459: {
463: 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);
464: if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
465: *compName = s->compNames[field][comp];
466: return(0);
467: }
469: /*@C
470: PetscSectionSetComponentName - Sets the name of a field component in the PetscSection
472: Not Collective
474: Input Parameters:
475: + s - the PetscSection
476: . field - the field number
477: . comp - the component number
478: - compName - the component name
480: Level: intermediate
482: .seealso: PetscSectionGetComponentName()
483: @*/
484: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
485: {
491: 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);
492: if ((comp < 0) || (comp >= s->numFieldComponents[field])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section field component %D should be in [%D, %D)", comp, 0, s->numFieldComponents[field]);
493: PetscFree(s->compNames[field][comp]);
494: PetscStrallocpy(compName, (char**) &s->compNames[field][comp]);
495: return(0);
496: }
498: /*@
499: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
501: Not collective
503: Input Parameters:
504: + s - the PetscSection
505: - field - the field number
507: Output Parameter:
508: . numComp - the number of field components
510: Level: intermediate
512: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
513: @*/
514: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
515: {
519: 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);
520: *numComp = s->numFieldComponents[field];
521: return(0);
522: }
524: /*@
525: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
527: Not collective
529: Input Parameters:
530: + s - the PetscSection
531: . field - the field number
532: - numComp - the number of field components
534: Level: intermediate
536: .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
537: @*/
538: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
539: {
541: PetscInt c;
542: char name[64];
546: 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);
547: if (s->compNames) {
548: for (c = 0; c < s->numFieldComponents[field]; ++c) {
549: PetscFree(s->compNames[field][c]);
550: }
551: PetscFree(s->compNames[field]);
552: }
554: s->numFieldComponents[field] = numComp;
555: if (numComp) {
556: PetscMalloc1(numComp, (char ***) &s->compNames[field]);
557: for (c = 0; c < numComp; ++c) {
558: PetscSNPrintf(name, 64, "%D", c);
559: PetscStrallocpy(name, (char **) &s->compNames[field][c]);
560: }
561: }
562: return(0);
563: }
565: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
566: {
570: if (!s->bc) {
571: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
572: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
573: }
574: return(0);
575: }
577: /*@
578: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
580: Not collective
582: Input Parameter:
583: . s - the PetscSection
585: Output Parameters:
586: + pStart - the first point
587: - pEnd - one past the last point
589: Level: intermediate
591: .seealso: PetscSectionSetChart(), PetscSectionCreate()
592: @*/
593: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
594: {
597: if (pStart) *pStart = s->pStart;
598: if (pEnd) *pEnd = s->pEnd;
599: return(0);
600: }
602: /*@
603: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
605: Not collective
607: Input Parameters:
608: + s - the PetscSection
609: . pStart - the first point
610: - pEnd - one past the last point
612: Level: intermediate
614: .seealso: PetscSectionGetChart(), PetscSectionCreate()
615: @*/
616: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
617: {
618: PetscInt f;
623: /* Cannot Reset() because it destroys field information */
624: s->setup = PETSC_FALSE;
625: PetscSectionDestroy(&s->bc);
626: PetscFree(s->bcIndices);
627: PetscFree2(s->atlasDof, s->atlasOff);
629: s->pStart = pStart;
630: s->pEnd = pEnd;
631: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
632: PetscArrayzero(s->atlasDof, pEnd - pStart);
633: for (f = 0; f < s->numFields; ++f) {
634: PetscSectionSetChart(s->field[f], pStart, pEnd);
635: }
636: return(0);
637: }
639: /*@
640: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
642: Not collective
644: Input Parameter:
645: . s - the PetscSection
647: Output Parameters:
648: . perm - The permutation as an IS
650: Level: intermediate
652: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
653: @*/
654: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
655: {
659: return(0);
660: }
662: /*@
663: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
665: Not collective
667: Input Parameters:
668: + s - the PetscSection
669: - perm - the permutation of points
671: Level: intermediate
673: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
674: @*/
675: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
676: {
682: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
683: if (s->perm != perm) {
684: ISDestroy(&s->perm);
685: if (perm) {
686: s->perm = perm;
687: PetscObjectReference((PetscObject) s->perm);
688: }
689: }
690: return(0);
691: }
693: /*@
694: PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major
696: Not collective
698: Input Parameter:
699: . s - the PetscSection
701: Output Parameter:
702: . pm - the flag for point major ordering
704: Level: intermediate
706: .seealso: PetscSectionSetPointMajor()
707: @*/
708: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
709: {
713: *pm = s->pointMajor;
714: return(0);
715: }
717: /*@
718: PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major
720: Not collective
722: Input Parameters:
723: + s - the PetscSection
724: - pm - the flag for point major ordering
726: Not collective
728: Level: intermediate
730: .seealso: PetscSectionGetPointMajor()
731: @*/
732: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
733: {
736: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
737: s->pointMajor = pm;
738: return(0);
739: }
741: /*@
742: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
744: Not collective
746: Input Parameters:
747: + s - the PetscSection
748: - point - the point
750: Output Parameter:
751: . numDof - the number of dof
753: Level: intermediate
755: .seealso: PetscSectionSetDof(), PetscSectionCreate()
756: @*/
757: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
758: {
762: #if defined(PETSC_USE_DEBUG)
763: 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);
764: #endif
765: *numDof = s->atlasDof[point - s->pStart];
766: return(0);
767: }
769: /*@
770: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
772: Not collective
774: Input Parameters:
775: + s - the PetscSection
776: . point - the point
777: - numDof - the number of dof
779: Level: intermediate
781: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
782: @*/
783: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
784: {
787: 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);
788: s->atlasDof[point - s->pStart] = numDof;
789: return(0);
790: }
792: /*@
793: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
795: Not collective
797: Input Parameters:
798: + s - the PetscSection
799: . point - the point
800: - numDof - the number of additional dof
802: Level: intermediate
804: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
805: @*/
806: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
807: {
810: #if defined(PETSC_USE_DEBUG)
811: 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);
812: #endif
813: s->atlasDof[point - s->pStart] += numDof;
814: return(0);
815: }
817: /*@
818: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
820: Not collective
822: Input Parameters:
823: + s - the PetscSection
824: . point - the point
825: - field - the field
827: Output Parameter:
828: . numDof - the number of dof
830: Level: intermediate
832: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
833: @*/
834: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
835: {
840: 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);
841: PetscSectionGetDof(s->field[field], point, numDof);
842: return(0);
843: }
845: /*@
846: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
848: Not collective
850: Input Parameters:
851: + s - the PetscSection
852: . point - the point
853: . field - the field
854: - numDof - the number of dof
856: Level: intermediate
858: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
859: @*/
860: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
861: {
866: 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);
867: PetscSectionSetDof(s->field[field], point, numDof);
868: return(0);
869: }
871: /*@
872: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
874: Not collective
876: Input Parameters:
877: + s - the PetscSection
878: . point - the point
879: . field - the field
880: - numDof - the number of dof
882: Level: intermediate
884: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
885: @*/
886: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
887: {
892: 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);
893: PetscSectionAddDof(s->field[field], point, numDof);
894: return(0);
895: }
897: /*@
898: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
900: Not collective
902: Input Parameters:
903: + s - the PetscSection
904: - point - the point
906: Output Parameter:
907: . numDof - the number of dof which are fixed by constraints
909: Level: intermediate
911: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
912: @*/
913: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
914: {
920: if (s->bc) {
921: PetscSectionGetDof(s->bc, point, numDof);
922: } else *numDof = 0;
923: return(0);
924: }
926: /*@
927: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
929: Not collective
931: Input Parameters:
932: + s - the PetscSection
933: . point - the point
934: - numDof - the number of dof which are fixed by constraints
936: Level: intermediate
938: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
939: @*/
940: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
941: {
946: if (numDof) {
947: PetscSectionCheckConstraints_Static(s);
948: PetscSectionSetDof(s->bc, point, numDof);
949: }
950: return(0);
951: }
953: /*@
954: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
956: Not collective
958: Input Parameters:
959: + s - the PetscSection
960: . point - the point
961: - numDof - the number of additional dof which are fixed by constraints
963: Level: intermediate
965: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
966: @*/
967: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
968: {
973: if (numDof) {
974: PetscSectionCheckConstraints_Static(s);
975: PetscSectionAddDof(s->bc, point, numDof);
976: }
977: return(0);
978: }
980: /*@
981: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
983: Not collective
985: Input Parameters:
986: + s - the PetscSection
987: . point - the point
988: - field - the field
990: Output Parameter:
991: . numDof - the number of dof which are fixed by constraints
993: Level: intermediate
995: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
996: @*/
997: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
998: {
1004: 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);
1005: PetscSectionGetConstraintDof(s->field[field], point, numDof);
1006: return(0);
1007: }
1009: /*@
1010: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1012: Not collective
1014: Input Parameters:
1015: + s - the PetscSection
1016: . point - the point
1017: . field - the field
1018: - numDof - the number of dof which are fixed by constraints
1020: Level: intermediate
1022: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1023: @*/
1024: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1025: {
1030: 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);
1031: PetscSectionSetConstraintDof(s->field[field], point, numDof);
1032: return(0);
1033: }
1035: /*@
1036: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1038: Not collective
1040: Input Parameters:
1041: + s - the PetscSection
1042: . point - the point
1043: . field - the field
1044: - numDof - the number of additional dof which are fixed by constraints
1046: Level: intermediate
1048: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
1049: @*/
1050: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1051: {
1056: 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);
1057: PetscSectionAddConstraintDof(s->field[field], point, numDof);
1058: return(0);
1059: }
1061: /*@
1062: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1064: Not collective
1066: Input Parameter:
1067: . s - the PetscSection
1069: Level: advanced
1071: .seealso: PetscSectionSetUp(), PetscSectionCreate()
1072: @*/
1073: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1074: {
1079: if (s->bc) {
1080: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
1082: PetscSectionSetUp(s->bc);
1083: PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1084: }
1085: return(0);
1086: }
1088: /*@
1089: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
1091: Not collective
1093: Input Parameter:
1094: . s - the PetscSection
1096: Level: intermediate
1098: .seealso: PetscSectionCreate()
1099: @*/
1100: PetscErrorCode PetscSectionSetUp(PetscSection s)
1101: {
1102: const PetscInt *pind = NULL;
1103: PetscInt offset = 0, foff, p, f;
1104: PetscErrorCode ierr;
1108: if (s->setup) return(0);
1109: s->setup = PETSC_TRUE;
1110: /* Set offsets and field offsets for all points */
1111: /* Assume that all fields have the same chart */
1112: if (s->perm) {ISGetIndices(s->perm, &pind);}
1113: if (s->pointMajor) {
1114: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1115: const PetscInt q = pind ? pind[p] : p;
1117: /* Set point offset */
1118: s->atlasOff[q] = offset;
1119: offset += s->atlasDof[q];
1120: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
1121: /* Set field offset */
1122: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1123: PetscSection sf = s->field[f];
1125: sf->atlasOff[q] = foff;
1126: foff += sf->atlasDof[q];
1127: }
1128: }
1129: } else {
1130: /* Set field offsets for all points */
1131: for (f = 0; f < s->numFields; ++f) {
1132: PetscSection sf = s->field[f];
1134: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1135: const PetscInt q = pind ? pind[p] : p;
1137: sf->atlasOff[q] = offset;
1138: offset += sf->atlasDof[q];
1139: }
1140: }
1141: /* Disable point offsets since these are unused */
1142: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1143: s->atlasOff[p] = -1;
1144: s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1145: }
1146: }
1147: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1148: /* Setup BC sections */
1149: PetscSectionSetUpBC(s);
1150: for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1151: return(0);
1152: }
1154: /*@
1155: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1157: Not collective
1159: Input Parameters:
1160: . s - the PetscSection
1162: Output Parameter:
1163: . maxDof - the maximum dof
1165: Level: intermediate
1167: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1168: @*/
1169: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1170: {
1174: *maxDof = s->maxDof;
1175: return(0);
1176: }
1178: /*@
1179: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
1181: Not collective
1183: Input Parameter:
1184: . s - the PetscSection
1186: Output Parameter:
1187: . size - the size of an array which can hold all the dofs
1189: Level: intermediate
1191: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1192: @*/
1193: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1194: {
1195: PetscInt p, n = 0;
1200: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1201: *size = n;
1202: return(0);
1203: }
1205: /*@
1206: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1208: Not collective
1210: Input Parameter:
1211: . s - the PetscSection
1213: Output Parameter:
1214: . size - the size of an array which can hold all unconstrained dofs
1216: Level: intermediate
1218: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1219: @*/
1220: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1221: {
1222: PetscInt p, n = 0;
1227: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1228: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1229: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1230: }
1231: *size = n;
1232: return(0);
1233: }
1235: /*@
1236: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1237: the local section and an SF describing the section point overlap.
1239: Input Parameters:
1240: + s - The PetscSection for the local field layout
1241: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1242: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1243: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1245: Output Parameter:
1246: . gsection - The PetscSection for the global field layout
1248: Note: This gives negative sizes and offsets to points not owned by this process
1250: Level: intermediate
1252: .seealso: PetscSectionCreate()
1253: @*/
1254: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1255: {
1256: PetscSection gs;
1257: const PetscInt *pind = NULL;
1258: PetscInt *recv = NULL, *neg = NULL;
1259: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1260: PetscErrorCode ierr;
1268: PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1269: PetscSectionGetChart(s, &pStart, &pEnd);
1270: PetscSectionSetChart(gs, pStart, pEnd);
1271: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1272: nlocal = nroots; /* The local/leaf space matches global/root space */
1273: /* Must allocate for all points visible to SF, which may be more than this section */
1274: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1275: PetscSFGetLeafRange(sf, NULL, &maxleaf);
1276: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1277: if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1278: PetscMalloc2(nroots,&neg,nlocal,&recv);
1279: PetscArrayzero(neg,nroots);
1280: }
1281: /* Mark all local points with negative dof */
1282: for (p = pStart; p < pEnd; ++p) {
1283: PetscSectionGetDof(s, p, &dof);
1284: PetscSectionSetDof(gs, p, dof);
1285: PetscSectionGetConstraintDof(s, p, &cdof);
1286: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1287: if (neg) neg[p] = -(dof+1);
1288: }
1289: PetscSectionSetUpBC(gs);
1290: if (gs->bcIndices) {PetscArraycpy(gs->bcIndices, s->bcIndices,gs->bc->atlasOff[gs->bc->pEnd-gs->bc->pStart-1] + gs->bc->atlasDof[gs->bc->pEnd-gs->bc->pStart-1]);}
1291: if (nroots >= 0) {
1292: PetscArrayzero(recv,nlocal);
1293: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1294: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1295: for (p = pStart; p < pEnd; ++p) {
1296: if (recv[p] < 0) {
1297: gs->atlasDof[p-pStart] = recv[p];
1298: PetscSectionGetDof(s, p, &dof);
1299: 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);
1300: }
1301: }
1302: }
1303: /* Calculate new sizes, get process offset, and calculate point offsets */
1304: if (s->perm) {ISGetIndices(s->perm, &pind);}
1305: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1306: const PetscInt q = pind ? pind[p] : p;
1308: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1309: gs->atlasOff[q] = off;
1310: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1311: }
1312: if (!localOffsets) {
1313: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1314: globalOff -= off;
1315: }
1316: for (p = pStart, off = 0; p < pEnd; ++p) {
1317: gs->atlasOff[p-pStart] += globalOff;
1318: if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1319: }
1320: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1321: /* Put in negative offsets for ghost points */
1322: if (nroots >= 0) {
1323: PetscArrayzero(recv,nlocal);
1324: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1325: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1326: for (p = pStart; p < pEnd; ++p) {
1327: if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1328: }
1329: }
1330: PetscFree2(neg,recv);
1331: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1332: *gsection = gs;
1333: return(0);
1334: }
1336: /*@
1337: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1338: the local section and an SF describing the section point overlap.
1340: Input Parameters:
1341: + s - The PetscSection for the local field layout
1342: . sf - The SF describing parallel layout of the section points
1343: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1344: . numExcludes - The number of exclusion ranges
1345: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1347: Output Parameter:
1348: . gsection - The PetscSection for the global field layout
1350: Note: This gives negative sizes and offsets to points not owned by this process
1352: Level: advanced
1354: .seealso: PetscSectionCreate()
1355: @*/
1356: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1357: {
1358: const PetscInt *pind = NULL;
1359: PetscInt *neg = NULL, *tmpOff = NULL;
1360: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1361: PetscErrorCode ierr;
1367: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1368: PetscSectionGetChart(s, &pStart, &pEnd);
1369: PetscSectionSetChart(*gsection, pStart, pEnd);
1370: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1371: if (nroots >= 0) {
1372: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1373: PetscCalloc1(nroots, &neg);
1374: if (nroots > pEnd-pStart) {
1375: PetscCalloc1(nroots, &tmpOff);
1376: } else {
1377: tmpOff = &(*gsection)->atlasDof[-pStart];
1378: }
1379: }
1380: /* Mark ghost points with negative dof */
1381: for (p = pStart; p < pEnd; ++p) {
1382: for (e = 0; e < numExcludes; ++e) {
1383: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1384: PetscSectionSetDof(*gsection, p, 0);
1385: break;
1386: }
1387: }
1388: if (e < numExcludes) continue;
1389: PetscSectionGetDof(s, p, &dof);
1390: PetscSectionSetDof(*gsection, p, dof);
1391: PetscSectionGetConstraintDof(s, p, &cdof);
1392: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1393: if (neg) neg[p] = -(dof+1);
1394: }
1395: PetscSectionSetUpBC(*gsection);
1396: if (nroots >= 0) {
1397: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1398: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1399: if (nroots > pEnd - pStart) {
1400: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1401: }
1402: }
1403: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1404: if (s->perm) {ISGetIndices(s->perm, &pind);}
1405: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1406: const PetscInt q = pind ? pind[p] : p;
1408: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1409: (*gsection)->atlasOff[q] = off;
1410: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1411: }
1412: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1413: globalOff -= off;
1414: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1415: (*gsection)->atlasOff[p] += globalOff;
1416: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1417: }
1418: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1419: /* Put in negative offsets for ghost points */
1420: if (nroots >= 0) {
1421: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1422: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1423: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1424: if (nroots > pEnd - pStart) {
1425: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1426: }
1427: }
1428: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1429: PetscFree(neg);
1430: return(0);
1431: }
1433: /*@
1434: PetscSectionGetPointLayout - Get the PetscLayout associated with the section points
1436: Collective on comm
1438: Input Parameters:
1439: + comm - The MPI_Comm
1440: - s - The PetscSection
1442: Output Parameter:
1443: . layout - The point layout for the section
1445: Note: This is usually caleld for the default global section.
1447: Level: advanced
1449: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1450: @*/
1451: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1452: {
1453: PetscInt pStart, pEnd, p, localSize = 0;
1457: PetscSectionGetChart(s, &pStart, &pEnd);
1458: for (p = pStart; p < pEnd; ++p) {
1459: PetscInt dof;
1461: PetscSectionGetDof(s, p, &dof);
1462: if (dof > 0) ++localSize;
1463: }
1464: PetscLayoutCreate(comm, layout);
1465: PetscLayoutSetLocalSize(*layout, localSize);
1466: PetscLayoutSetBlockSize(*layout, 1);
1467: PetscLayoutSetUp(*layout);
1468: return(0);
1469: }
1471: /*@
1472: PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.
1474: Collective on comm
1476: Input Parameters:
1477: + comm - The MPI_Comm
1478: - s - The PetscSection
1480: Output Parameter:
1481: . layout - The dof layout for the section
1483: Note: This is usually called for the default global section.
1485: Level: advanced
1487: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1488: @*/
1489: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1490: {
1491: PetscInt pStart, pEnd, p, localSize = 0;
1497: PetscSectionGetChart(s, &pStart, &pEnd);
1498: for (p = pStart; p < pEnd; ++p) {
1499: PetscInt dof,cdof;
1501: PetscSectionGetDof(s, p, &dof);
1502: PetscSectionGetConstraintDof(s, p, &cdof);
1503: if (dof-cdof > 0) localSize += dof-cdof;
1504: }
1505: PetscLayoutCreate(comm, layout);
1506: PetscLayoutSetLocalSize(*layout, localSize);
1507: PetscLayoutSetBlockSize(*layout, 1);
1508: PetscLayoutSetUp(*layout);
1509: return(0);
1510: }
1512: /*@
1513: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1515: Not collective
1517: Input Parameters:
1518: + s - the PetscSection
1519: - point - the point
1521: Output Parameter:
1522: . offset - the offset
1524: Level: intermediate
1526: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1527: @*/
1528: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1529: {
1533: #if defined(PETSC_USE_DEBUG)
1534: 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);
1535: #endif
1536: *offset = s->atlasOff[point - s->pStart];
1537: return(0);
1538: }
1540: /*@
1541: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1543: Not collective
1545: Input Parameters:
1546: + s - the PetscSection
1547: . point - the point
1548: - offset - the offset
1550: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1552: Level: intermediate
1554: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1555: @*/
1556: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1557: {
1560: 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);
1561: s->atlasOff[point - s->pStart] = offset;
1562: return(0);
1563: }
1565: /*@
1566: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1568: Not collective
1570: Input Parameters:
1571: + s - the PetscSection
1572: . point - the point
1573: - field - the field
1575: Output Parameter:
1576: . offset - the offset
1578: Level: intermediate
1580: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1581: @*/
1582: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1583: {
1589: 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);
1590: PetscSectionGetOffset(s->field[field], point, offset);
1591: return(0);
1592: }
1594: /*@
1595: PetscSectionSetFieldOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1597: Not collective
1599: Input Parameters:
1600: + s - the PetscSection
1601: . point - the point
1602: . field - the field
1603: - offset - the offset
1605: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1607: Level: intermediate
1609: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1610: @*/
1611: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1612: {
1617: 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);
1618: PetscSectionSetOffset(s->field[field], point, offset);
1619: return(0);
1620: }
1622: /*@
1623: PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.
1625: Not collective
1627: Input Parameters:
1628: + s - the PetscSection
1629: . point - the point
1630: - field - the field
1632: Output Parameter:
1633: . offset - the offset
1635: Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1636: this point, what number is the first dof with this field.
1638: Level: advanced
1640: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1641: @*/
1642: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1643: {
1644: PetscInt off, foff;
1650: 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);
1651: PetscSectionGetOffset(s, point, &off);
1652: PetscSectionGetOffset(s->field[field], point, &foff);
1653: *offset = foff - off;
1654: return(0);
1655: }
1657: /*@
1658: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1660: Not collective
1662: Input Parameter:
1663: . s - the PetscSection
1665: Output Parameters:
1666: + start - the minimum offset
1667: - end - one more than the maximum offset
1669: Level: intermediate
1671: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1672: @*/
1673: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1674: {
1675: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1680: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1681: PetscSectionGetChart(s, &pStart, &pEnd);
1682: for (p = 0; p < pEnd-pStart; ++p) {
1683: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1685: if (off >= 0) {
1686: os = PetscMin(os, off);
1687: oe = PetscMax(oe, off+dof);
1688: }
1689: }
1690: if (start) *start = os;
1691: if (end) *end = oe;
1692: return(0);
1693: }
1695: /*@
1696: PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields
1698: Collective on s
1700: Input Parameter:
1701: + s - the PetscSection
1702: . len - the number of subfields
1703: - fields - the subfield numbers
1705: Output Parameter:
1706: . subs - the subsection
1708: Note: The section offsets now refer to a new, smaller vector.
1710: Level: advanced
1712: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1713: @*/
1714: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1715: {
1716: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
1720: if (!len) return(0);
1724: PetscSectionGetNumFields(s, &nF);
1725: if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1726: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1727: PetscSectionSetNumFields(*subs, len);
1728: for (f = 0; f < len; ++f) {
1729: const char *name = NULL;
1730: PetscInt numComp = 0;
1732: PetscSectionGetFieldName(s, fields[f], &name);
1733: PetscSectionSetFieldName(*subs, f, name);
1734: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1735: PetscSectionSetFieldComponents(*subs, f, numComp);
1736: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
1737: PetscSectionGetComponentName(s, fields[f], c, &name);
1738: PetscSectionSetComponentName(*subs, f, c, name);
1739: }
1740: }
1741: PetscSectionGetChart(s, &pStart, &pEnd);
1742: PetscSectionSetChart(*subs, pStart, pEnd);
1743: for (p = pStart; p < pEnd; ++p) {
1744: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1746: for (f = 0; f < len; ++f) {
1747: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1748: PetscSectionSetFieldDof(*subs, p, f, fdof);
1749: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1750: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1751: dof += fdof;
1752: cdof += cfdof;
1753: }
1754: PetscSectionSetDof(*subs, p, dof);
1755: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1756: maxCdof = PetscMax(cdof, maxCdof);
1757: }
1758: PetscSectionSetUp(*subs);
1759: if (maxCdof) {
1760: PetscInt *indices;
1762: PetscMalloc1(maxCdof, &indices);
1763: for (p = pStart; p < pEnd; ++p) {
1764: PetscInt cdof;
1766: PetscSectionGetConstraintDof(*subs, p, &cdof);
1767: if (cdof) {
1768: const PetscInt *oldIndices = NULL;
1769: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1771: for (f = 0; f < len; ++f) {
1772: PetscInt oldFoff = 0, oldf;
1774: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1775: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1776: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1777: /* This can be sped up if we assume sorted fields */
1778: for (oldf = 0; oldf < fields[f]; ++oldf) {
1779: PetscInt oldfdof = 0;
1780: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1781: oldFoff += oldfdof;
1782: }
1783: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1784: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1785: numConst += cfdof;
1786: fOff += fdof;
1787: }
1788: PetscSectionSetConstraintIndices(*subs, p, indices);
1789: }
1790: }
1791: PetscFree(indices);
1792: }
1793: return(0);
1794: }
1796: /*@
1797: PetscSectionCreateSupersection - Create a new, larger section composed of the input sections
1799: Collective on s
1801: Input Parameters:
1802: + s - the input sections
1803: - len - the number of input sections
1805: Output Parameter:
1806: . supers - the supersection
1808: Note: The section offsets now refer to a new, larger vector.
1810: Level: advanced
1812: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1813: @*/
1814: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1815: {
1816: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1820: if (!len) return(0);
1821: for (i = 0; i < len; ++i) {
1822: PetscInt nf, pStarti, pEndi;
1824: PetscSectionGetNumFields(s[i], &nf);
1825: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1826: pStart = PetscMin(pStart, pStarti);
1827: pEnd = PetscMax(pEnd, pEndi);
1828: Nf += nf;
1829: }
1830: PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1831: PetscSectionSetNumFields(*supers, Nf);
1832: for (i = 0, f = 0; i < len; ++i) {
1833: PetscInt nf, fi, ci;
1835: PetscSectionGetNumFields(s[i], &nf);
1836: for (fi = 0; fi < nf; ++fi, ++f) {
1837: const char *name = NULL;
1838: PetscInt numComp = 0;
1840: PetscSectionGetFieldName(s[i], fi, &name);
1841: PetscSectionSetFieldName(*supers, f, name);
1842: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1843: PetscSectionSetFieldComponents(*supers, f, numComp);
1844: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
1845: PetscSectionGetComponentName(s[i], fi, ci, &name);
1846: PetscSectionSetComponentName(*supers, f, ci, name);
1847: }
1848: }
1849: }
1850: PetscSectionSetChart(*supers, pStart, pEnd);
1851: for (p = pStart; p < pEnd; ++p) {
1852: PetscInt dof = 0, cdof = 0;
1854: for (i = 0, f = 0; i < len; ++i) {
1855: PetscInt nf, fi, pStarti, pEndi;
1856: PetscInt fdof = 0, cfdof = 0;
1858: PetscSectionGetNumFields(s[i], &nf);
1859: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1860: if ((p < pStarti) || (p >= pEndi)) continue;
1861: for (fi = 0; fi < nf; ++fi, ++f) {
1862: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1863: PetscSectionAddFieldDof(*supers, p, f, fdof);
1864: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1865: if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1866: dof += fdof;
1867: cdof += cfdof;
1868: }
1869: }
1870: PetscSectionSetDof(*supers, p, dof);
1871: if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1872: maxCdof = PetscMax(cdof, maxCdof);
1873: }
1874: PetscSectionSetUp(*supers);
1875: if (maxCdof) {
1876: PetscInt *indices;
1878: PetscMalloc1(maxCdof, &indices);
1879: for (p = pStart; p < pEnd; ++p) {
1880: PetscInt cdof;
1882: PetscSectionGetConstraintDof(*supers, p, &cdof);
1883: if (cdof) {
1884: PetscInt dof, numConst = 0, fOff = 0;
1886: for (i = 0, f = 0; i < len; ++i) {
1887: const PetscInt *oldIndices = NULL;
1888: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1890: PetscSectionGetNumFields(s[i], &nf);
1891: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1892: if ((p < pStarti) || (p >= pEndi)) continue;
1893: for (fi = 0; fi < nf; ++fi, ++f) {
1894: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1895: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1896: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1897: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1898: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1899: numConst += cfdof;
1900: }
1901: PetscSectionGetDof(s[i], p, &dof);
1902: fOff += dof;
1903: }
1904: PetscSectionSetConstraintIndices(*supers, p, indices);
1905: }
1906: }
1907: PetscFree(indices);
1908: }
1909: return(0);
1910: }
1912: /*@
1913: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
1915: Collective on s
1917: Input Parameters:
1918: + s - the PetscSection
1919: - subpointMap - a sorted list of points in the original mesh which are in the submesh
1921: Output Parameter:
1922: . subs - the subsection
1924: Note: The section offsets now refer to a new, smaller vector.
1926: Level: advanced
1928: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1929: @*/
1930: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1931: {
1932: const PetscInt *points = NULL, *indices = NULL;
1933: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, subp;
1940: PetscSectionGetNumFields(s, &numFields);
1941: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1942: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1943: for (f = 0; f < numFields; ++f) {
1944: const char *name = NULL;
1945: PetscInt numComp = 0;
1947: PetscSectionGetFieldName(s, f, &name);
1948: PetscSectionSetFieldName(*subs, f, name);
1949: PetscSectionGetFieldComponents(s, f, &numComp);
1950: PetscSectionSetFieldComponents(*subs, f, numComp);
1951: for (c = 0; c < s->numFieldComponents[f]; ++c) {
1952: PetscSectionGetComponentName(s, f, c, &name);
1953: PetscSectionSetComponentName(*subs, f, c, name);
1954: }
1955: }
1956: /* For right now, we do not try to squeeze the subchart */
1957: if (subpointMap) {
1958: ISGetSize(subpointMap, &numSubpoints);
1959: ISGetIndices(subpointMap, &points);
1960: }
1961: PetscSectionGetChart(s, &pStart, &pEnd);
1962: PetscSectionSetChart(*subs, 0, numSubpoints);
1963: for (p = pStart; p < pEnd; ++p) {
1964: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1966: PetscFindInt(p, numSubpoints, points, &subp);
1967: if (subp < 0) continue;
1968: for (f = 0; f < numFields; ++f) {
1969: PetscSectionGetFieldDof(s, p, f, &fdof);
1970: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1971: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1972: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1973: }
1974: PetscSectionGetDof(s, p, &dof);
1975: PetscSectionSetDof(*subs, subp, dof);
1976: PetscSectionGetConstraintDof(s, p, &cdof);
1977: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1978: }
1979: PetscSectionSetUp(*subs);
1980: /* Change offsets to original offsets */
1981: for (p = pStart; p < pEnd; ++p) {
1982: PetscInt off, foff = 0;
1984: PetscFindInt(p, numSubpoints, points, &subp);
1985: if (subp < 0) continue;
1986: for (f = 0; f < numFields; ++f) {
1987: PetscSectionGetFieldOffset(s, p, f, &foff);
1988: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1989: }
1990: PetscSectionGetOffset(s, p, &off);
1991: PetscSectionSetOffset(*subs, subp, off);
1992: }
1993: /* Copy constraint indices */
1994: for (subp = 0; subp < numSubpoints; ++subp) {
1995: PetscInt cdof;
1997: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1998: if (cdof) {
1999: for (f = 0; f < numFields; ++f) {
2000: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
2001: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
2002: }
2003: PetscSectionGetConstraintIndices(s, points[subp], &indices);
2004: PetscSectionSetConstraintIndices(*subs, subp, indices);
2005: }
2006: }
2007: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
2008: return(0);
2009: }
2011: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2012: {
2013: PetscInt p;
2014: PetscMPIInt rank;
2018: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
2019: PetscViewerASCIIPushSynchronized(viewer);
2020: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
2021: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2022: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
2023: PetscInt b;
2025: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2026: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
2027: PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
2028: }
2029: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
2030: } else {
2031: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
2032: }
2033: }
2034: PetscViewerFlush(viewer);
2035: PetscViewerASCIIPopSynchronized(viewer);
2036: if (s->sym) {
2037: PetscViewerASCIIPushTab(viewer);
2038: PetscSectionSymView(s->sym,viewer);
2039: PetscViewerASCIIPopTab(viewer);
2040: }
2041: return(0);
2042: }
2044: /*@C
2045: PetscSectionViewFromOptions - View from Options
2047: Collective on PetscSection
2049: Input Parameters:
2050: + A - the PetscSection object to view
2051: . obj - Optional object
2052: - name - command line option
2054: Level: intermediate
2055: .seealso: PetscSection, PetscSectionView, PetscObjectViewFromOptions(), PetscSectionCreate()
2056: @*/
2057: PetscErrorCode PetscSectionViewFromOptions(PetscSection A,PetscObject obj,const char name[])
2058: {
2063: PetscObjectViewFromOptions((PetscObject)A,obj,name);
2064: return(0);
2065: }
2067: /*@C
2068: PetscSectionView - Views a PetscSection
2070: Collective on PetscSection
2072: Input Parameters:
2073: + s - the PetscSection object to view
2074: - v - the viewer
2076: Level: beginner
2078: .seealso PetscSectionCreate(), PetscSectionDestroy()
2079: @*/
2080: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2081: {
2082: PetscBool isascii;
2083: PetscInt f;
2088: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
2090: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
2091: if (isascii) {
2092: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
2093: if (s->numFields) {
2094: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
2095: for (f = 0; f < s->numFields; ++f) {
2096: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
2097: PetscSectionView_ASCII(s->field[f], viewer);
2098: }
2099: } else {
2100: PetscSectionView_ASCII(s, viewer);
2101: }
2102: }
2103: return(0);
2104: }
2106: /*@
2107: PetscSectionReset - Frees all section data.
2109: Not collective
2111: Input Parameters:
2112: . s - the PetscSection
2114: Level: beginner
2116: .seealso: PetscSection, PetscSectionCreate()
2117: @*/
2118: PetscErrorCode PetscSectionReset(PetscSection s)
2119: {
2120: PetscInt f, c;
2125: for (f = 0; f < s->numFields; ++f) {
2126: PetscSectionDestroy(&s->field[f]);
2127: PetscFree(s->fieldNames[f]);
2128: for (c = 0; c < s->numFieldComponents[f]; ++c)
2129: PetscFree(s->compNames[f][c]);
2130: PetscFree(s->compNames[f]);
2131: }
2132: PetscFree(s->numFieldComponents);
2133: PetscFree(s->fieldNames);
2134: PetscFree(s->compNames);
2135: PetscFree(s->field);
2136: PetscSectionDestroy(&s->bc);
2137: PetscFree(s->bcIndices);
2138: PetscFree2(s->atlasDof, s->atlasOff);
2139: PetscSectionDestroy(&s->clSection);
2140: ISDestroy(&s->clPoints);
2141: ISDestroy(&s->perm);
2142: PetscFree(s->clPerm);
2143: PetscFree(s->clInvPerm);
2144: PetscSectionSymDestroy(&s->sym);
2145: PetscSectionDestroy(&s->clSection);
2146: ISDestroy(&s->clPoints);
2148: s->pStart = -1;
2149: s->pEnd = -1;
2150: s->maxDof = 0;
2151: s->setup = PETSC_FALSE;
2152: s->numFields = 0;
2153: s->clObj = NULL;
2154: return(0);
2155: }
2157: /*@
2158: PetscSectionDestroy - Frees a section object and frees its range if that exists.
2160: Not collective
2162: Input Parameters:
2163: . s - the PetscSection
2165: Level: beginner
2167: .seealso: PetscSection, PetscSectionCreate()
2168: @*/
2169: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2170: {
2174: if (!*s) return(0);
2176: if (--((PetscObject)(*s))->refct > 0) {
2177: *s = NULL;
2178: return(0);
2179: }
2180: PetscSectionReset(*s);
2181: PetscHeaderDestroy(s);
2182: return(0);
2183: }
2185: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2186: {
2187: const PetscInt p = point - s->pStart;
2191: *values = &baseArray[s->atlasOff[p]];
2192: return(0);
2193: }
2195: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2196: {
2197: PetscInt *array;
2198: const PetscInt p = point - s->pStart;
2199: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2200: PetscInt cDim = 0;
2205: PetscSectionGetConstraintDof(s, p, &cDim);
2206: array = &baseArray[s->atlasOff[p]];
2207: if (!cDim) {
2208: if (orientation >= 0) {
2209: const PetscInt dim = s->atlasDof[p];
2210: PetscInt i;
2212: if (mode == INSERT_VALUES) {
2213: for (i = 0; i < dim; ++i) array[i] = values[i];
2214: } else {
2215: for (i = 0; i < dim; ++i) array[i] += values[i];
2216: }
2217: } else {
2218: PetscInt offset = 0;
2219: PetscInt j = -1, field, i;
2221: for (field = 0; field < s->numFields; ++field) {
2222: const PetscInt dim = s->field[field]->atlasDof[p];
2224: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2225: offset += dim;
2226: }
2227: }
2228: } else {
2229: if (orientation >= 0) {
2230: const PetscInt dim = s->atlasDof[p];
2231: PetscInt cInd = 0, i;
2232: const PetscInt *cDof;
2234: PetscSectionGetConstraintIndices(s, point, &cDof);
2235: if (mode == INSERT_VALUES) {
2236: for (i = 0; i < dim; ++i) {
2237: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2238: array[i] = values[i];
2239: }
2240: } else {
2241: for (i = 0; i < dim; ++i) {
2242: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2243: array[i] += values[i];
2244: }
2245: }
2246: } else {
2247: const PetscInt *cDof;
2248: PetscInt offset = 0;
2249: PetscInt cOffset = 0;
2250: PetscInt j = 0, field;
2252: PetscSectionGetConstraintIndices(s, point, &cDof);
2253: for (field = 0; field < s->numFields; ++field) {
2254: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2255: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2256: const PetscInt sDim = dim - tDim;
2257: PetscInt cInd = 0, i,k;
2259: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2260: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2261: array[j] = values[k];
2262: }
2263: offset += dim;
2264: cOffset += dim - tDim;
2265: }
2266: }
2267: }
2268: return(0);
2269: }
2271: /*@C
2272: PetscSectionHasConstraints - Determine whether a section has constrained dofs
2274: Not collective
2276: Input Parameter:
2277: . s - The PetscSection
2279: Output Parameter:
2280: . hasConstraints - flag indicating that the section has constrained dofs
2282: Level: intermediate
2284: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2285: @*/
2286: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2287: {
2291: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2292: return(0);
2293: }
2295: /*@C
2296: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2298: Not collective
2300: Input Parameters:
2301: + s - The PetscSection
2302: - point - The point
2304: Output Parameter:
2305: . indices - The constrained dofs
2307: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
2309: Level: intermediate
2311: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2312: @*/
2313: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2314: {
2319: if (s->bc) {
2320: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2321: } else *indices = NULL;
2322: return(0);
2323: }
2325: /*@C
2326: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2328: Not collective
2330: Input Parameters:
2331: + s - The PetscSection
2332: . point - The point
2333: - indices - The constrained dofs
2335: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
2337: Level: intermediate
2339: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2340: @*/
2341: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2342: {
2347: if (s->bc) {
2348: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2349: }
2350: return(0);
2351: }
2353: /*@C
2354: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2356: Not collective
2358: Input Parameters:
2359: + s - The PetscSection
2360: . field - The field number
2361: - point - The point
2363: Output Parameter:
2364: . indices - The constrained dofs sorted in ascending order
2366: Notes:
2367: The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by PetscSectionGetConstraintDof().
2369: Fortran Note:
2370: In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()
2372: Level: intermediate
2374: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2375: @*/
2376: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2377: {
2382: 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);
2383: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2384: return(0);
2385: }
2387: /*@C
2388: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2390: Not collective
2392: Input Parameters:
2393: + s - The PetscSection
2394: . point - The point
2395: . field - The field number
2396: - indices - The constrained dofs
2398: Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()
2400: Level: intermediate
2402: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2403: @*/
2404: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2405: {
2410: 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);
2411: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2412: return(0);
2413: }
2415: /*@
2416: PetscSectionPermute - Reorder the section according to the input point permutation
2418: Collective on PetscSection
2420: Input Parameter:
2421: + section - The PetscSection object
2422: - perm - The point permutation, old point p becomes new point perm[p]
2424: Output Parameter:
2425: . sectionNew - The permuted PetscSection
2427: Level: intermediate
2429: .seealso: MatPermute()
2430: @*/
2431: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2432: {
2433: PetscSection s = section, sNew;
2434: const PetscInt *perm;
2435: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2436: PetscErrorCode ierr;
2442: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2443: PetscSectionGetNumFields(s, &numFields);
2444: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2445: for (f = 0; f < numFields; ++f) {
2446: const char *name;
2447: PetscInt numComp;
2449: PetscSectionGetFieldName(s, f, &name);
2450: PetscSectionSetFieldName(sNew, f, name);
2451: PetscSectionGetFieldComponents(s, f, &numComp);
2452: PetscSectionSetFieldComponents(sNew, f, numComp);
2453: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2454: PetscSectionGetComponentName(s, f, c, &name);
2455: PetscSectionSetComponentName(sNew, f, c, name);
2456: }
2457: }
2458: ISGetLocalSize(permutation, &numPoints);
2459: ISGetIndices(permutation, &perm);
2460: PetscSectionGetChart(s, &pStart, &pEnd);
2461: PetscSectionSetChart(sNew, pStart, pEnd);
2462: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2463: for (p = pStart; p < pEnd; ++p) {
2464: PetscInt dof, cdof;
2466: PetscSectionGetDof(s, p, &dof);
2467: PetscSectionSetDof(sNew, perm[p], dof);
2468: PetscSectionGetConstraintDof(s, p, &cdof);
2469: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2470: for (f = 0; f < numFields; ++f) {
2471: PetscSectionGetFieldDof(s, p, f, &dof);
2472: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2473: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2474: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2475: }
2476: }
2477: PetscSectionSetUp(sNew);
2478: for (p = pStart; p < pEnd; ++p) {
2479: const PetscInt *cind;
2480: PetscInt cdof;
2482: PetscSectionGetConstraintDof(s, p, &cdof);
2483: if (cdof) {
2484: PetscSectionGetConstraintIndices(s, p, &cind);
2485: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2486: }
2487: for (f = 0; f < numFields; ++f) {
2488: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2489: if (cdof) {
2490: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2491: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2492: }
2493: }
2494: }
2495: ISRestoreIndices(permutation, &perm);
2496: *sectionNew = sNew;
2497: return(0);
2498: }
2500: /* TODO: the next three functions should be moved to sf/utils */
2501: #include <petsc/private/sfimpl.h>
2503: /*@C
2504: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
2506: Collective on sf
2508: Input Parameters:
2509: + sf - The SF
2510: - rootSection - Section defined on root space
2512: Output Parameters:
2513: + remoteOffsets - root offsets in leaf storage, or NULL
2514: - leafSection - Section defined on the leaf space
2516: Level: advanced
2518: .seealso: PetscSFCreate()
2519: @*/
2520: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2521: {
2522: PetscSF embedSF;
2523: const PetscInt *indices;
2524: IS selected;
2525: PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
2526: PetscBool *sub, hasc;
2530: PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
2531: PetscSectionGetNumFields(rootSection, &numFields);
2532: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2533: PetscMalloc1(numFields+2, &sub);
2534: sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
2535: for (f = 0; f < numFields; ++f) {
2536: PetscSectionSym sym;
2537: const char *name = NULL;
2538: PetscInt numComp = 0;
2540: sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
2541: PetscSectionGetFieldComponents(rootSection, f, &numComp);
2542: PetscSectionGetFieldName(rootSection, f, &name);
2543: PetscSectionGetFieldSym(rootSection, f, &sym);
2544: PetscSectionSetFieldComponents(leafSection, f, numComp);
2545: PetscSectionSetFieldName(leafSection, f, name);
2546: PetscSectionSetFieldSym(leafSection, f, sym);
2547: for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
2548: PetscSectionGetComponentName(rootSection, f, c, &name);
2549: PetscSectionSetComponentName(leafSection, f, c, name);
2550: }
2551: }
2552: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2553: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2554: rpEnd = PetscMin(rpEnd,nroots);
2555: rpEnd = PetscMax(rpStart,rpEnd);
2556: /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
2557: sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
2558: MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
2559: if (sub[0]) {
2560: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2561: ISGetIndices(selected, &indices);
2562: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2563: ISRestoreIndices(selected, &indices);
2564: ISDestroy(&selected);
2565: } else {
2566: PetscObjectReference((PetscObject)sf);
2567: embedSF = sf;
2568: }
2569: PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
2570: lpEnd++;
2572: PetscSectionSetChart(leafSection, lpStart, lpEnd);
2574: /* Constrained dof section */
2575: hasc = sub[1];
2576: for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);
2578: /* Could fuse these at the cost of copies and extra allocation */
2579: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2580: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2581: if (sub[1]) {
2582: PetscSectionCheckConstraints_Static(rootSection);
2583: PetscSectionCheckConstraints_Static(leafSection);
2584: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2585: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2586: }
2587: for (f = 0; f < numFields; ++f) {
2588: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2589: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2590: if (sub[2+f]) {
2591: PetscSectionCheckConstraints_Static(rootSection->field[f]);
2592: PetscSectionCheckConstraints_Static(leafSection->field[f]);
2593: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2594: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2595: }
2596: }
2597: if (remoteOffsets) {
2598: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2599: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2600: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2601: }
2602: PetscSectionSetUp(leafSection);
2603: if (hasc) { /* need to communicate bcIndices */
2604: PetscSF bcSF;
2605: PetscInt *rOffBc;
2607: PetscMalloc1(lpEnd - lpStart, &rOffBc);
2608: if (sub[1]) {
2609: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2610: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2611: PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
2612: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2613: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2614: PetscSFDestroy(&bcSF);
2615: }
2616: for (f = 0; f < numFields; ++f) {
2617: if (sub[2+f]) {
2618: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2619: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2620: PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
2621: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2622: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2623: PetscSFDestroy(&bcSF);
2624: }
2625: }
2626: PetscFree(rOffBc);
2627: }
2628: PetscSFDestroy(&embedSF);
2629: PetscFree(sub);
2630: PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
2631: return(0);
2632: }
2634: /*@C
2635: PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
2637: Collective on sf
2639: Input Parameters:
2640: + sf - The SF
2641: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
2642: - leafSection - Data layout of local points for incoming data (this is layout for SF leaves)
2644: Output Parameter:
2645: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2647: Level: developer
2649: .seealso: PetscSFCreate()
2650: @*/
2651: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2652: {
2653: PetscSF embedSF;
2654: const PetscInt *indices;
2655: IS selected;
2656: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2657: PetscErrorCode ierr;
2660: *remoteOffsets = NULL;
2661: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2662: if (numRoots < 0) return(0);
2663: PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
2664: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2665: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2666: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2667: ISGetIndices(selected, &indices);
2668: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2669: ISRestoreIndices(selected, &indices);
2670: ISDestroy(&selected);
2671: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2672: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2673: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2674: PetscSFDestroy(&embedSF);
2675: PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
2676: return(0);
2677: }
2679: /*@C
2680: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2682: Collective on sf
2684: Input Parameters:
2685: + sf - The SF
2686: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2687: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2688: - leafSection - Data layout of local points for incoming data (this is the distributed section)
2690: Output Parameters:
2691: - sectionSF - The new SF
2693: Note: Either rootSection or remoteOffsets can be specified
2695: Level: advanced
2697: .seealso: PetscSFCreate()
2698: @*/
2699: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2700: {
2701: MPI_Comm comm;
2702: const PetscInt *localPoints;
2703: const PetscSFNode *remotePoints;
2704: PetscInt lpStart, lpEnd;
2705: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2706: PetscInt *localIndices;
2707: PetscSFNode *remoteIndices;
2708: PetscInt i, ind;
2709: PetscErrorCode ierr;
2717: PetscObjectGetComm((PetscObject)sf,&comm);
2718: PetscSFCreate(comm, sectionSF);
2719: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2720: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2721: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2722: if (numRoots < 0) return(0);
2723: PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
2724: for (i = 0; i < numPoints; ++i) {
2725: PetscInt localPoint = localPoints ? localPoints[i] : i;
2726: PetscInt dof;
2728: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2729: PetscSectionGetDof(leafSection, localPoint, &dof);
2730: numIndices += dof;
2731: }
2732: }
2733: PetscMalloc1(numIndices, &localIndices);
2734: PetscMalloc1(numIndices, &remoteIndices);
2735: /* Create new index graph */
2736: for (i = 0, ind = 0; i < numPoints; ++i) {
2737: PetscInt localPoint = localPoints ? localPoints[i] : i;
2738: PetscInt rank = remotePoints[i].rank;
2740: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2741: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2742: PetscInt loff, dof, d;
2744: PetscSectionGetOffset(leafSection, localPoint, &loff);
2745: PetscSectionGetDof(leafSection, localPoint, &dof);
2746: for (d = 0; d < dof; ++d, ++ind) {
2747: localIndices[ind] = loff+d;
2748: remoteIndices[ind].rank = rank;
2749: remoteIndices[ind].index = remoteOffset+d;
2750: }
2751: }
2752: }
2753: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2754: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2755: PetscSFSetUp(*sectionSF);
2756: PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
2757: return(0);
2758: }
2760: /*@
2761: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2763: Collective on section
2765: Input Parameters:
2766: + section - The PetscSection
2767: . obj - A PetscObject which serves as the key for this index
2768: . clSection - Section giving the size of the closure of each point
2769: - clPoints - IS giving the points in each closure
2771: Note: We compress out closure points with no dofs in this section
2773: Level: advanced
2775: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2776: @*/
2777: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2778: {
2785: if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2786: section->clObj = obj;
2787: PetscObjectReference((PetscObject)clSection);
2788: PetscObjectReference((PetscObject)clPoints);
2789: PetscSectionDestroy(§ion->clSection);
2790: ISDestroy(§ion->clPoints);
2791: section->clSection = clSection;
2792: section->clPoints = clPoints;
2793: return(0);
2794: }
2796: /*@
2797: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2799: Collective on section
2801: Input Parameters:
2802: + section - The PetscSection
2803: - obj - A PetscObject which serves as the key for this index
2805: Output Parameters:
2806: + clSection - Section giving the size of the closure of each point
2807: - clPoints - IS giving the points in each closure
2809: Note: We compress out closure points with no dofs in this section
2811: Level: advanced
2813: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2814: @*/
2815: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2816: {
2818: if (section->clObj == obj) {
2819: if (clSection) *clSection = section->clSection;
2820: if (clPoints) *clPoints = section->clPoints;
2821: } else {
2822: if (clSection) *clSection = NULL;
2823: if (clPoints) *clPoints = NULL;
2824: }
2825: return(0);
2826: }
2828: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2829: {
2830: PetscInt i;
2834: if (section->clObj != obj) {
2835: PetscSectionDestroy(§ion->clSection);
2836: ISDestroy(§ion->clPoints);
2837: }
2838: section->clObj = obj;
2839: PetscFree(section->clPerm);
2840: PetscFree(section->clInvPerm);
2841: section->clSize = clSize;
2842: if (mode == PETSC_COPY_VALUES) {
2843: PetscMalloc1(clSize, §ion->clPerm);
2844: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2845: PetscArraycpy(section->clPerm, clPerm, clSize);
2846: } else if (mode == PETSC_OWN_POINTER) {
2847: section->clPerm = clPerm;
2848: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2849: PetscMalloc1(clSize, §ion->clInvPerm);
2850: for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2851: return(0);
2852: }
2854: /*@
2855: PetscSectionSetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2857: Not Collective
2859: Input Parameters:
2860: + section - The PetscSection
2861: . obj - A PetscObject which serves as the key for this index
2862: - perm - Permutation of the cell dof closure
2864: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2865: other points (like faces).
2867: Level: intermediate
2869: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2870: @*/
2871: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2872: {
2873: const PetscInt *clPerm = NULL;
2874: PetscInt clSize = 0;
2875: PetscErrorCode ierr;
2878: if (perm) {
2879: ISGetLocalSize(perm, &clSize);
2880: ISGetIndices(perm, &clPerm);
2881: }
2882: PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2883: if (perm) {ISRestoreIndices(perm, &clPerm);}
2884: return(0);
2885: }
2887: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2888: {
2890: if (section->clObj == obj) {
2891: if (size) *size = section->clSize;
2892: if (perm) *perm = section->clPerm;
2893: } else {
2894: if (size) *size = 0;
2895: if (perm) *perm = NULL;
2896: }
2897: return(0);
2898: }
2900: /*@
2901: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2903: Not collective
2905: Input Parameters:
2906: + section - The PetscSection
2907: - obj - A PetscObject which serves as the key for this index
2909: Output Parameter:
2910: . perm - The dof closure permutation
2912: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2913: other points (like faces).
2915: The user must destroy the IS that is returned.
2917: Level: intermediate
2919: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2920: @*/
2921: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2922: {
2923: const PetscInt *clPerm;
2924: PetscInt clSize;
2925: PetscErrorCode ierr;
2928: PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2929: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2930: return(0);
2931: }
2933: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2934: {
2936: if (section->clObj == obj) {
2937: if (size) *size = section->clSize;
2938: if (perm) *perm = section->clInvPerm;
2939: } else {
2940: if (size) *size = 0;
2941: if (perm) *perm = NULL;
2942: }
2943: return(0);
2944: }
2946: /*@
2947: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2949: Not collective
2951: Input Parameters:
2952: + section - The PetscSection
2953: - obj - A PetscObject which serves as the key for this index
2955: Output Parameters:
2956: + size - The dof closure size
2957: - perm - The dof closure permutation
2959: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2960: other points (like faces).
2962: The user must destroy the IS that is returned.
2964: Level: intermediate
2966: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2967: @*/
2968: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2969: {
2970: const PetscInt *clPerm;
2971: PetscInt clSize;
2972: PetscErrorCode ierr;
2975: PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2976: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2977: return(0);
2978: }
2980: /*@
2981: PetscSectionGetField - Get the subsection associated with a single field
2983: Input Parameters:
2984: + s - The PetscSection
2985: - field - The field number
2987: Output Parameter:
2988: . subs - The subsection for the given field
2990: Level: intermediate
2992: .seealso: PetscSectionSetNumFields()
2993: @*/
2994: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2995: {
2999: 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);
3000: *subs = s->field[field];
3001: return(0);
3002: }
3004: PetscClassId PETSC_SECTION_SYM_CLASSID;
3005: PetscFunctionList PetscSectionSymList = NULL;
3007: /*@
3008: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
3010: Collective
3012: Input Parameter:
3013: . comm - the MPI communicator
3015: Output Parameter:
3016: . sym - pointer to the new set of symmetries
3018: Level: developer
3020: .seealso: PetscSectionSym, PetscSectionSymDestroy()
3021: @*/
3022: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3023: {
3028: ISInitializePackage();
3029: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
3030: return(0);
3031: }
3033: /*@C
3034: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
3036: Collective on PetscSectionSym
3038: Input Parameters:
3039: + sym - The section symmetry object
3040: - method - The name of the section symmetry type
3042: Level: developer
3044: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
3045: @*/
3046: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3047: {
3048: PetscErrorCode (*r)(PetscSectionSym);
3049: PetscBool match;
3054: PetscObjectTypeCompare((PetscObject) sym, method, &match);
3055: if (match) return(0);
3057: PetscFunctionListFind(PetscSectionSymList,method,&r);
3058: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3059: if (sym->ops->destroy) {
3060: (*sym->ops->destroy)(sym);
3061: sym->ops->destroy = NULL;
3062: }
3063: (*r)(sym);
3064: PetscObjectChangeTypeName((PetscObject)sym,method);
3065: return(0);
3066: }
3069: /*@C
3070: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
3072: Not Collective
3074: Input Parameter:
3075: . sym - The section symmetry
3077: Output Parameter:
3078: . type - The index set type name
3080: Level: developer
3082: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
3083: @*/
3084: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3085: {
3089: *type = ((PetscObject)sym)->type_name;
3090: return(0);
3091: }
3093: /*@C
3094: PetscSectionSymRegister - Adds a new section symmetry implementation
3096: Not Collective
3098: Input Parameters:
3099: + name - The name of a new user-defined creation routine
3100: - create_func - The creation routine itself
3102: Notes:
3103: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
3105: Level: developer
3107: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
3108: @*/
3109: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3110: {
3114: ISInitializePackage();
3115: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
3116: return(0);
3117: }
3119: /*@
3120: PetscSectionSymDestroy - Destroys a section symmetry.
3122: Collective on PetscSectionSym
3124: Input Parameters:
3125: . sym - the section symmetry
3127: Level: developer
3129: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
3130: @*/
3131: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3132: {
3133: SymWorkLink link,next;
3137: if (!*sym) return(0);
3139: if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
3140: if ((*sym)->ops->destroy) {
3141: (*(*sym)->ops->destroy)(*sym);
3142: }
3143: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3144: for (link=(*sym)->workin; link; link=next) {
3145: next = link->next;
3146: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
3147: PetscFree(link);
3148: }
3149: (*sym)->workin = NULL;
3150: PetscHeaderDestroy(sym);
3151: return(0);
3152: }
3154: /*@C
3155: PetscSectionSymView - Displays a section symmetry
3157: Collective on PetscSectionSym
3159: Input Parameters:
3160: + sym - the index set
3161: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
3163: Level: developer
3165: .seealso: PetscViewerASCIIOpen()
3166: @*/
3167: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3168: {
3173: if (!viewer) {
3174: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
3175: }
3178: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3179: if (sym->ops->view) {
3180: (*sym->ops->view)(sym,viewer);
3181: }
3182: return(0);
3183: }
3185: /*@
3186: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3188: Collective on PetscSection
3190: Input Parameters:
3191: + section - the section describing data layout
3192: - sym - the symmetry describing the affect of orientation on the access of the data
3194: Level: developer
3196: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3197: @*/
3198: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3199: {
3204: PetscSectionSymDestroy(&(section->sym));
3205: if (sym) {
3208: PetscObjectReference((PetscObject) sym);
3209: }
3210: section->sym = sym;
3211: return(0);
3212: }
3214: /*@
3215: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3217: Not collective
3219: Input Parameters:
3220: . section - the section describing data layout
3222: Output Parameters:
3223: . sym - the symmetry describing the affect of orientation on the access of the data
3225: Level: developer
3227: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3228: @*/
3229: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3230: {
3233: *sym = section->sym;
3234: return(0);
3235: }
3237: /*@
3238: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3240: Collective on PetscSection
3242: Input Parameters:
3243: + section - the section describing data layout
3244: . field - the field number
3245: - sym - the symmetry describing the affect of orientation on the access of the data
3247: Level: developer
3249: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3250: @*/
3251: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3252: {
3257: 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);
3258: PetscSectionSetSym(section->field[field],sym);
3259: return(0);
3260: }
3262: /*@
3263: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3265: Collective on PetscSection
3267: Input Parameters:
3268: + section - the section describing data layout
3269: - field - the field number
3271: Output Parameters:
3272: . sym - the symmetry describing the affect of orientation on the access of the data
3274: Level: developer
3276: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3277: @*/
3278: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3279: {
3282: 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);
3283: *sym = section->field[field]->sym;
3284: return(0);
3285: }
3287: /*@C
3288: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
3290: Not collective
3292: Input Parameters:
3293: + section - the section
3294: . numPoints - the number of points
3295: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3296: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3297: context, see DMPlexGetConeOrientation()).
3299: Output Parameter:
3300: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3301: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3302: identity).
3304: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3305: .vb
3306: const PetscInt **perms;
3307: const PetscScalar **rots;
3308: PetscInt lOffset;
3310: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3311: for (i = 0, lOffset = 0; i < numPoints; i++) {
3312: PetscInt point = points[2*i], dof, sOffset;
3313: const PetscInt *perm = perms ? perms[i] : NULL;
3314: const PetscScalar *rot = rots ? rots[i] : NULL;
3316: PetscSectionGetDof(section,point,&dof);
3317: PetscSectionGetOffset(section,point,&sOffset);
3319: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3320: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3321: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3322: lOffset += dof;
3323: }
3324: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3325: .ve
3327: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3328: .vb
3329: const PetscInt **perms;
3330: const PetscScalar **rots;
3331: PetscInt lOffset;
3333: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3334: for (i = 0, lOffset = 0; i < numPoints; i++) {
3335: PetscInt point = points[2*i], dof, sOffset;
3336: const PetscInt *perm = perms ? perms[i] : NULL;
3337: const PetscScalar *rot = rots ? rots[i] : NULL;
3339: PetscSectionGetDof(section,point,&dof);
3340: PetscSectionGetOffset(section,point,&sOff);
3342: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3343: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3344: offset += dof;
3345: }
3346: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3347: .ve
3349: Level: developer
3351: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3352: @*/
3353: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3354: {
3355: PetscSectionSym sym;
3356: PetscErrorCode ierr;
3361: if (perms) *perms = NULL;
3362: if (rots) *rots = NULL;
3363: sym = section->sym;
3364: if (sym && (perms || rots)) {
3365: SymWorkLink link;
3367: if (sym->workin) {
3368: link = sym->workin;
3369: sym->workin = sym->workin->next;
3370: } else {
3371: PetscNewLog(sym,&link);
3372: }
3373: if (numPoints > link->numPoints) {
3374: PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3375: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3376: link->numPoints = numPoints;
3377: }
3378: link->next = sym->workout;
3379: sym->workout = link;
3380: PetscArrayzero((PetscInt**)link->perms,numPoints);
3381: PetscArrayzero((PetscInt**)link->rots,numPoints);
3382: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3383: if (perms) *perms = link->perms;
3384: if (rots) *rots = link->rots;
3385: }
3386: return(0);
3387: }
3389: /*@C
3390: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
3392: Not collective
3394: Input Parameters:
3395: + section - the section
3396: . numPoints - the number of points
3397: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3398: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3399: context, see DMPlexGetConeOrientation()).
3401: Output Parameter:
3402: + perms - The permutations for the given orientations: set to NULL at conclusion
3403: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3405: Level: developer
3407: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3408: @*/
3409: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3410: {
3411: PetscSectionSym sym;
3415: sym = section->sym;
3416: if (sym && (perms || rots)) {
3417: SymWorkLink *p,link;
3419: for (p=&sym->workout; (link=*p); p=&link->next) {
3420: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3421: *p = link->next;
3422: link->next = sym->workin;
3423: sym->workin = link;
3424: if (perms) *perms = NULL;
3425: if (rots) *rots = NULL;
3426: return(0);
3427: }
3428: }
3429: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3430: }
3431: return(0);
3432: }
3434: /*@C
3435: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
3437: Not collective
3439: Input Parameters:
3440: + section - the section
3441: . field - the field of the section
3442: . numPoints - the number of points
3443: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3444: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3445: context, see DMPlexGetConeOrientation()).
3447: Output Parameter:
3448: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3449: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3450: identity).
3452: Level: developer
3454: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3455: @*/
3456: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3457: {
3462: 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);
3463: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3464: return(0);
3465: }
3467: /*@C
3468: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
3470: Not collective
3472: Input Parameters:
3473: + section - the section
3474: . field - the field number
3475: . numPoints - the number of points
3476: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3477: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3478: context, see DMPlexGetConeOrientation()).
3480: Output Parameter:
3481: + perms - The permutations for the given orientations: set to NULL at conclusion
3482: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3484: Level: developer
3486: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3487: @*/
3488: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3489: {
3494: if (field > section->numFields) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"field %D greater than number of fields (%D) in section",field,section->numFields);
3495: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3496: return(0);
3497: }
3499: /*@
3500: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3502: Not collective
3504: Input Parameter:
3505: . s - the global PetscSection
3507: Output Parameters:
3508: . flg - the flag
3510: Level: developer
3512: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3513: @*/
3514: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3515: {
3518: *flg = s->useFieldOff;
3519: return(0);
3520: }
3522: /*@
3523: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3525: Not collective
3527: Input Parameters:
3528: + s - the global PetscSection
3529: - flg - the flag
3531: Level: developer
3533: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3534: @*/
3535: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3536: {
3539: s->useFieldOff = flg;
3540: return(0);
3541: }
3543: #define PetscSectionExpandPoints_Loop(TYPE) \
3544: { \
3545: PetscInt i, n, o0, o1, size; \
3546: TYPE *a0 = (TYPE*)origArray, *a1; \
3547: PetscSectionGetStorageSize(s, &size); \
3548: PetscMalloc1(size, &a1); \
3549: for (i=0; i<npoints; i++) { \
3550: PetscSectionGetOffset(origSection, points_[i], &o0); \
3551: PetscSectionGetOffset(s, i, &o1); \
3552: PetscSectionGetDof(s, i, &n); \
3553: PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3554: } \
3555: *newArray = (void*)a1; \
3556: }
3558: /*@
3559: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3561: Not collective
3563: Input Parameters:
3564: + origSection - the PetscSection describing the layout of the array
3565: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3566: . origArray - the array; its size must be equal to the storage size of origSection
3567: - points - IS with points to extract; its indices must lie in the chart of origSection
3569: Output Parameters:
3570: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3571: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3573: Level: developer
3575: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3576: @*/
3577: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3578: {
3579: PetscSection s;
3580: const PetscInt *points_;
3581: PetscInt i, n, npoints, pStart, pEnd;
3582: PetscMPIInt unitsize;
3583: PetscErrorCode ierr;
3591: MPI_Type_size(dataType, &unitsize);
3592: ISGetLocalSize(points, &npoints);
3593: ISGetIndices(points, &points_);
3594: PetscSectionGetChart(origSection, &pStart, &pEnd);
3595: PetscSectionCreate(PETSC_COMM_SELF, &s);
3596: PetscSectionSetChart(s, 0, npoints);
3597: for (i=0; i<npoints; i++) {
3598: 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);
3599: PetscSectionGetDof(origSection, points_[i], &n);
3600: PetscSectionSetDof(s, i, n);
3601: }
3602: PetscSectionSetUp(s);
3603: if (newArray) {
3604: if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);}
3605: else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);}
3606: else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);}
3607: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3608: }
3609: if (newSection) {
3610: *newSection = s;
3611: } else {
3612: PetscSectionDestroy(&s);
3613: }
3614: return(0);
3615: }