Actual source code: section.c
petsc-3.12.5 2020-03-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)->clObj = NULL;
63: (*s)->clSection = NULL;
64: (*s)->clPoints = NULL;
65: (*s)->clSize = 0;
66: (*s)->clPerm = NULL;
67: (*s)->clInvPerm = NULL;
68: return(0);
69: }
71: /*@
72: PetscSectionCopy - Creates a shallow (if possible) copy of the PetscSection
74: Collective
76: Input Parameter:
77: . section - the PetscSection
79: Output Parameter:
80: . newSection - the copy
82: Level: intermediate
84: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
85: @*/
86: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
87: {
88: PetscSectionSym sym;
89: IS perm;
90: PetscInt numFields, f, pStart, pEnd, p;
91: PetscErrorCode ierr;
96: PetscSectionGetNumFields(section, &numFields);
97: if (numFields) {PetscSectionSetNumFields(newSection, numFields);}
98: for (f = 0; f < numFields; ++f) {
99: const char *name = NULL;
100: PetscInt numComp = 0;
102: PetscSectionGetFieldName(section, f, &name);
103: PetscSectionSetFieldName(newSection, f, name);
104: PetscSectionGetFieldComponents(section, f, &numComp);
105: PetscSectionSetFieldComponents(newSection, f, numComp);
106: PetscSectionGetFieldSym(section, f, &sym);
107: PetscSectionSetFieldSym(newSection, f, sym);
108: }
109: PetscSectionGetChart(section, &pStart, &pEnd);
110: PetscSectionSetChart(newSection, pStart, pEnd);
111: PetscSectionGetPermutation(section, &perm);
112: PetscSectionSetPermutation(newSection, perm);
113: PetscSectionGetSym(section, &sym);
114: PetscSectionSetSym(newSection, sym);
115: for (p = pStart; p < pEnd; ++p) {
116: PetscInt dof, cdof, fcdof = 0;
118: PetscSectionGetDof(section, p, &dof);
119: PetscSectionSetDof(newSection, p, dof);
120: PetscSectionGetConstraintDof(section, p, &cdof);
121: if (cdof) {PetscSectionSetConstraintDof(newSection, p, cdof);}
122: for (f = 0; f < numFields; ++f) {
123: PetscSectionGetFieldDof(section, p, f, &dof);
124: PetscSectionSetFieldDof(newSection, p, f, dof);
125: if (cdof) {
126: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
127: if (fcdof) {PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof);}
128: }
129: }
130: }
131: PetscSectionSetUp(newSection);
132: for (p = pStart; p < pEnd; ++p) {
133: PetscInt off, cdof, fcdof = 0;
134: const PetscInt *cInd;
136: /* Must set offsets in case they do not agree with the prefix sums */
137: PetscSectionGetOffset(section, p, &off);
138: PetscSectionSetOffset(newSection, p, off);
139: PetscSectionGetConstraintDof(section, p, &cdof);
140: if (cdof) {
141: PetscSectionGetConstraintIndices(section, p, &cInd);
142: PetscSectionSetConstraintIndices(newSection, p, cInd);
143: for (f = 0; f < numFields; ++f) {
144: PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
145: if (fcdof) {
146: PetscSectionGetFieldConstraintIndices(section, p, f, &cInd);
147: PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd);
148: }
149: }
150: }
151: }
152: return(0);
153: }
155: /*@
156: PetscSectionClone - Creates a shallow (if possible) copy of the PetscSection
158: Collective
160: Input Parameter:
161: . section - the PetscSection
163: Output Parameter:
164: . newSection - the copy
166: Level: beginner
168: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
169: @*/
170: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
171: {
177: PetscSectionCreate(PetscObjectComm((PetscObject) section), newSection);
178: PetscSectionCopy(section, *newSection);
179: return(0);
180: }
182: /*@
183: PetscSectionSetFromOptions - sets parameters in a PetscSection from the options database
185: Collective on PetscSection
187: Input Parameter:
188: . section - the PetscSection
190: Options Database:
191: . -petscsection_point_major the dof order
193: Level: intermediate
195: .seealso: PetscSection, PetscSectionCreate(), PetscSectionDestroy()
196: @*/
197: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
198: {
203: PetscObjectOptionsBegin((PetscObject) s);
204: PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL);
205: /* process any options handlers added with PetscObjectAddOptionsHandler() */
206: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) s);
207: PetscOptionsEnd();
208: PetscObjectViewFromOptions((PetscObject) s, NULL, "-petscsection_view");
209: return(0);
210: }
212: /*@
213: PetscSectionCompare - Compares two sections
215: Collective on PetscSection
217: Input Parameterss:
218: + s1 - the first PetscSection
219: - s2 - the second PetscSection
221: Output Parameter:
222: . congruent - PETSC_TRUE if the two sections are congruent, PETSC_FALSE otherwise
224: Level: intermediate
226: Notes:
227: Field names are disregarded.
229: .seealso: PetscSection, PetscSectionCreate(), PetscSectionCopy(), PetscSectionClone()
230: @*/
231: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
232: {
233: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
234: const PetscInt *idx1, *idx2;
235: IS perm1, perm2;
236: PetscBool flg;
237: PetscMPIInt mflg;
244: flg = PETSC_FALSE;
246: MPI_Comm_compare(PetscObjectComm((PetscObject)s1),PetscObjectComm((PetscObject)s2),&mflg);
247: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
248: *congruent = PETSC_FALSE;
249: return(0);
250: }
252: PetscSectionGetChart(s1, &pStart, &pEnd);
253: PetscSectionGetChart(s2, &n1, &n2);
254: if (pStart != n1 || pEnd != n2) goto not_congruent;
256: PetscSectionGetPermutation(s1, &perm1);
257: PetscSectionGetPermutation(s2, &perm2);
258: if (perm1 && perm2) {
259: ISEqual(perm1, perm2, congruent);
260: if (!(*congruent)) goto not_congruent;
261: } else if (perm1 != perm2) goto not_congruent;
263: for (p = pStart; p < pEnd; ++p) {
264: PetscSectionGetOffset(s1, p, &n1);
265: PetscSectionGetOffset(s2, p, &n2);
266: if (n1 != n2) goto not_congruent;
268: PetscSectionGetDof(s1, p, &n1);
269: PetscSectionGetDof(s2, p, &n2);
270: if (n1 != n2) goto not_congruent;
272: PetscSectionGetConstraintDof(s1, p, &ncdof);
273: PetscSectionGetConstraintDof(s2, p, &n2);
274: if (ncdof != n2) goto not_congruent;
276: PetscSectionGetConstraintIndices(s1, p, &idx1);
277: PetscSectionGetConstraintIndices(s2, p, &idx2);
278: PetscArraycmp(idx1, idx2, ncdof, congruent);
279: if (!(*congruent)) goto not_congruent;
280: }
282: PetscSectionGetNumFields(s1, &nfields);
283: PetscSectionGetNumFields(s2, &n2);
284: if (nfields != n2) goto not_congruent;
286: for (f = 0; f < nfields; ++f) {
287: PetscSectionGetFieldComponents(s1, f, &n1);
288: PetscSectionGetFieldComponents(s2, f, &n2);
289: if (n1 != n2) goto not_congruent;
291: for (p = pStart; p < pEnd; ++p) {
292: PetscSectionGetFieldOffset(s1, p, f, &n1);
293: PetscSectionGetFieldOffset(s2, p, f, &n2);
294: if (n1 != n2) goto not_congruent;
296: PetscSectionGetFieldDof(s1, p, f, &n1);
297: PetscSectionGetFieldDof(s2, p, f, &n2);
298: if (n1 != n2) goto not_congruent;
300: PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof);
301: PetscSectionGetFieldConstraintDof(s2, p, f, &n2);
302: if (nfcdof != n2) goto not_congruent;
304: PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1);
305: PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2);
306: PetscArraycmp(idx1, idx2, nfcdof, congruent);
307: if (!(*congruent)) goto not_congruent;
308: }
309: }
311: flg = PETSC_TRUE;
312: not_congruent:
313: MPIU_Allreduce(&flg,congruent,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)s1));
314: return(0);
315: }
317: /*@
318: PetscSectionGetNumFields - Returns the number of fields, or 0 if no fields were defined.
320: Not collective
322: Input Parameter:
323: . s - the PetscSection
325: Output Parameter:
326: . numFields - the number of fields defined, or 0 if none were defined
328: Level: intermediate
330: .seealso: PetscSectionSetNumFields()
331: @*/
332: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
333: {
337: *numFields = s->numFields;
338: return(0);
339: }
341: /*@
342: PetscSectionSetNumFields - Sets the number of fields.
344: Not collective
346: Input Parameters:
347: + s - the PetscSection
348: - numFields - the number of fields
350: Level: intermediate
352: .seealso: PetscSectionGetNumFields()
353: @*/
354: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
355: {
356: PetscInt f;
361: if (numFields <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %D must be positive", numFields);
362: PetscSectionReset(s);
364: s->numFields = numFields;
365: PetscMalloc1(s->numFields, &s->numFieldComponents);
366: PetscMalloc1(s->numFields, &s->fieldNames);
367: PetscMalloc1(s->numFields, &s->field);
368: for (f = 0; f < s->numFields; ++f) {
369: char name[64];
371: s->numFieldComponents[f] = 1;
373: PetscSectionCreate(PetscObjectComm((PetscObject) s), &s->field[f]);
374: PetscSNPrintf(name, 64, "Field_%D", f);
375: PetscStrallocpy(name, (char **) &s->fieldNames[f]);
376: }
377: return(0);
378: }
380: /*@C
381: PetscSectionGetFieldName - Returns the name of a field in the PetscSection
383: Not Collective
385: Input Parameters:
386: + s - the PetscSection
387: - field - the field number
389: Output Parameter:
390: . fieldName - the field name
392: Level: intermediate
394: .seealso: PetscSectionSetFieldName()
395: @*/
396: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
397: {
401: 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);
402: *fieldName = s->fieldNames[field];
403: return(0);
404: }
406: /*@C
407: PetscSectionSetFieldName - Sets the name of a field in the PetscSection
409: Not Collective
411: Input Parameters:
412: + s - the PetscSection
413: . field - the field number
414: - fieldName - the field name
416: Level: intermediate
418: .seealso: PetscSectionGetFieldName()
419: @*/
420: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
421: {
427: 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);
428: PetscFree(s->fieldNames[field]);
429: PetscStrallocpy(fieldName, (char**) &s->fieldNames[field]);
430: return(0);
431: }
433: /*@
434: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
436: Not collective
438: Input Parameters:
439: + s - the PetscSection
440: - field - the field number
442: Output Parameter:
443: . numComp - the number of field components
445: Level: intermediate
447: .seealso: PetscSectionSetFieldComponents(), PetscSectionGetNumFields()
448: @*/
449: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
450: {
454: 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);
455: *numComp = s->numFieldComponents[field];
456: return(0);
457: }
459: /*@
460: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
462: Not collective
464: Input Parameters:
465: + s - the PetscSection
466: . field - the field number
467: - numComp - the number of field components
469: Level: intermediate
471: .seealso: PetscSectionGetFieldComponents(), PetscSectionGetNumFields()
472: @*/
473: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
474: {
477: 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);
478: s->numFieldComponents[field] = numComp;
479: return(0);
480: }
482: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
483: {
487: if (!s->bc) {
488: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
489: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
490: }
491: return(0);
492: }
494: /*@
495: PetscSectionGetChart - Returns the range [pStart, pEnd) in which points lie.
497: Not collective
499: Input Parameter:
500: . s - the PetscSection
502: Output Parameters:
503: + pStart - the first point
504: - pEnd - one past the last point
506: Level: intermediate
508: .seealso: PetscSectionSetChart(), PetscSectionCreate()
509: @*/
510: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
511: {
514: if (pStart) *pStart = s->pStart;
515: if (pEnd) *pEnd = s->pEnd;
516: return(0);
517: }
519: /*@
520: PetscSectionSetChart - Sets the range [pStart, pEnd) in which points lie.
522: Not collective
524: Input Parameters:
525: + s - the PetscSection
526: . pStart - the first point
527: - pEnd - one past the last point
529: Level: intermediate
531: .seealso: PetscSectionGetChart(), PetscSectionCreate()
532: @*/
533: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
534: {
535: PetscInt f;
540: /* Cannot Reset() because it destroys field information */
541: s->setup = PETSC_FALSE;
542: PetscSectionDestroy(&s->bc);
543: PetscFree(s->bcIndices);
544: PetscFree2(s->atlasDof, s->atlasOff);
546: s->pStart = pStart;
547: s->pEnd = pEnd;
548: PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff);
549: PetscArrayzero(s->atlasDof, pEnd - pStart);
550: for (f = 0; f < s->numFields; ++f) {
551: PetscSectionSetChart(s->field[f], pStart, pEnd);
552: }
553: return(0);
554: }
556: /*@
557: PetscSectionGetPermutation - Returns the permutation of [0, pEnd-pStart) or NULL
559: Not collective
561: Input Parameter:
562: . s - the PetscSection
564: Output Parameters:
565: . perm - The permutation as an IS
567: Level: intermediate
569: .seealso: PetscSectionSetPermutation(), PetscSectionCreate()
570: @*/
571: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
572: {
576: return(0);
577: }
579: /*@
580: PetscSectionSetPermutation - Sets the permutation for [0, pEnd-pStart)
582: Not collective
584: Input Parameters:
585: + s - the PetscSection
586: - perm - the permutation of points
588: Level: intermediate
590: .seealso: PetscSectionGetPermutation(), PetscSectionCreate()
591: @*/
592: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
593: {
599: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
600: if (s->perm != perm) {
601: ISDestroy(&s->perm);
602: if (perm) {
603: s->perm = perm;
604: PetscObjectReference((PetscObject) s->perm);
605: }
606: }
607: return(0);
608: }
610: /*@
611: PetscSectionGetPointMajor - Returns the flag for dof ordering, true if it is point major, otherwise field major
613: Not collective
615: Input Parameter:
616: . s - the PetscSection
618: Output Parameter:
619: . pm - the flag for point major ordering
621: Level: intermediate
623: .seealso: PetscSectionSetPointMajor()
624: @*/
625: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
626: {
630: *pm = s->pointMajor;
631: return(0);
632: }
634: /*@
635: PetscSectionSetPointMajor - Sets the flag for dof ordering, true if it is point major, otherwise field major
637: Not collective
639: Input Parameters:
640: + s - the PetscSection
641: - pm - the flag for point major ordering
643: Not collective
645: Level: intermediate
647: .seealso: PetscSectionGetPointMajor()
648: @*/
649: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
650: {
653: if (s->setup) SETERRQ(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
654: s->pointMajor = pm;
655: return(0);
656: }
658: /*@
659: PetscSectionGetDof - Return the number of degrees of freedom associated with a given point.
661: Not collective
663: Input Parameters:
664: + s - the PetscSection
665: - point - the point
667: Output Parameter:
668: . numDof - the number of dof
670: Level: intermediate
672: .seealso: PetscSectionSetDof(), PetscSectionCreate()
673: @*/
674: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
675: {
679: #if defined(PETSC_USE_DEBUG)
680: 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);
681: #endif
682: *numDof = s->atlasDof[point - s->pStart];
683: return(0);
684: }
686: /*@
687: PetscSectionSetDof - Sets the number of degrees of freedom associated with a given point.
689: Not collective
691: Input Parameters:
692: + s - the PetscSection
693: . point - the point
694: - numDof - the number of dof
696: Level: intermediate
698: .seealso: PetscSectionGetDof(), PetscSectionAddDof(), PetscSectionCreate()
699: @*/
700: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
701: {
704: 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);
705: s->atlasDof[point - s->pStart] = numDof;
706: return(0);
707: }
709: /*@
710: PetscSectionAddDof - Adds to the number of degrees of freedom associated with a given point.
712: Not collective
714: Input Parameters:
715: + s - the PetscSection
716: . point - the point
717: - numDof - the number of additional dof
719: Level: intermediate
721: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
722: @*/
723: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
724: {
727: #if defined(PETSC_USE_DEBUG)
728: 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);
729: #endif
730: s->atlasDof[point - s->pStart] += numDof;
731: return(0);
732: }
734: /*@
735: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
737: Not collective
739: Input Parameters:
740: + s - the PetscSection
741: . point - the point
742: - field - the field
744: Output Parameter:
745: . numDof - the number of dof
747: Level: intermediate
749: .seealso: PetscSectionSetFieldDof(), PetscSectionCreate()
750: @*/
751: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
752: {
757: 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);
758: PetscSectionGetDof(s->field[field], point, numDof);
759: return(0);
760: }
762: /*@
763: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
765: Not collective
767: Input Parameters:
768: + s - the PetscSection
769: . point - the point
770: . field - the field
771: - numDof - the number of dof
773: Level: intermediate
775: .seealso: PetscSectionGetFieldDof(), PetscSectionCreate()
776: @*/
777: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
778: {
783: 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);
784: PetscSectionSetDof(s->field[field], point, numDof);
785: return(0);
786: }
788: /*@
789: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
791: Not collective
793: Input Parameters:
794: + s - the PetscSection
795: . point - the point
796: . field - the field
797: - numDof - the number of dof
799: Level: intermediate
801: .seealso: PetscSectionSetFieldDof(), PetscSectionGetFieldDof(), PetscSectionCreate()
802: @*/
803: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
804: {
809: 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);
810: PetscSectionAddDof(s->field[field], point, numDof);
811: return(0);
812: }
814: /*@
815: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
817: Not collective
819: Input Parameters:
820: + s - the PetscSection
821: - point - the point
823: Output Parameter:
824: . numDof - the number of dof which are fixed by constraints
826: Level: intermediate
828: .seealso: PetscSectionGetDof(), PetscSectionSetConstraintDof(), PetscSectionCreate()
829: @*/
830: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
831: {
837: if (s->bc) {
838: PetscSectionGetDof(s->bc, point, numDof);
839: } else *numDof = 0;
840: return(0);
841: }
843: /*@
844: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
846: Not collective
848: Input Parameters:
849: + s - the PetscSection
850: . point - the point
851: - numDof - the number of dof which are fixed by constraints
853: Level: intermediate
855: .seealso: PetscSectionSetDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
856: @*/
857: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
858: {
863: if (numDof) {
864: PetscSectionCheckConstraints_Static(s);
865: PetscSectionSetDof(s->bc, point, numDof);
866: }
867: return(0);
868: }
870: /*@
871: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
873: Not collective
875: Input Parameters:
876: + s - the PetscSection
877: . point - the point
878: - numDof - the number of additional dof which are fixed by constraints
880: Level: intermediate
882: .seealso: PetscSectionAddDof(), PetscSectionGetConstraintDof(), PetscSectionCreate()
883: @*/
884: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
885: {
890: if (numDof) {
891: PetscSectionCheckConstraints_Static(s);
892: PetscSectionAddDof(s->bc, point, numDof);
893: }
894: return(0);
895: }
897: /*@
898: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
900: Not collective
902: Input Parameters:
903: + s - the PetscSection
904: . point - the point
905: - field - the field
907: Output Parameter:
908: . numDof - the number of dof which are fixed by constraints
910: Level: intermediate
912: .seealso: PetscSectionGetDof(), PetscSectionSetFieldConstraintDof(), PetscSectionCreate()
913: @*/
914: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
915: {
921: 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);
922: PetscSectionGetConstraintDof(s->field[field], point, numDof);
923: return(0);
924: }
926: /*@
927: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
929: Not collective
931: Input Parameters:
932: + s - the PetscSection
933: . point - the point
934: . field - the field
935: - numDof - the number of dof which are fixed by constraints
937: Level: intermediate
939: .seealso: PetscSectionSetDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
940: @*/
941: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
942: {
947: 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);
948: PetscSectionSetConstraintDof(s->field[field], point, numDof);
949: return(0);
950: }
952: /*@
953: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
955: Not collective
957: Input Parameters:
958: + s - the PetscSection
959: . point - the point
960: . field - the field
961: - numDof - the number of additional dof which are fixed by constraints
963: Level: intermediate
965: .seealso: PetscSectionAddDof(), PetscSectionGetFieldConstraintDof(), PetscSectionCreate()
966: @*/
967: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
968: {
973: 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);
974: PetscSectionAddConstraintDof(s->field[field], point, numDof);
975: return(0);
976: }
978: /*@
979: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
981: Not collective
983: Input Parameter:
984: . s - the PetscSection
986: Level: advanced
988: .seealso: PetscSectionSetUp(), PetscSectionCreate()
989: @*/
990: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
991: {
996: if (s->bc) {
997: const PetscInt last = (s->bc->pEnd-s->bc->pStart) - 1;
999: PetscSectionSetUp(s->bc);
1000: PetscMalloc1(last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0, &s->bcIndices);
1001: }
1002: return(0);
1003: }
1005: /*@
1006: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point.
1008: Not collective
1010: Input Parameter:
1011: . s - the PetscSection
1013: Level: intermediate
1015: .seealso: PetscSectionCreate()
1016: @*/
1017: PetscErrorCode PetscSectionSetUp(PetscSection s)
1018: {
1019: const PetscInt *pind = NULL;
1020: PetscInt offset = 0, foff, p, f;
1021: PetscErrorCode ierr;
1025: if (s->setup) return(0);
1026: s->setup = PETSC_TRUE;
1027: /* Set offsets and field offsets for all points */
1028: /* Assume that all fields have the same chart */
1029: if (s->perm) {ISGetIndices(s->perm, &pind);}
1030: if (s->pointMajor) {
1031: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1032: const PetscInt q = pind ? pind[p] : p;
1034: /* Set point offset */
1035: s->atlasOff[q] = offset;
1036: offset += s->atlasDof[q];
1037: s->maxDof = PetscMax(s->maxDof, s->atlasDof[q]);
1038: /* Set field offset */
1039: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1040: PetscSection sf = s->field[f];
1042: sf->atlasOff[q] = foff;
1043: foff += sf->atlasDof[q];
1044: }
1045: }
1046: } else {
1047: /* Set field offsets for all points */
1048: for (f = 0; f < s->numFields; ++f) {
1049: PetscSection sf = s->field[f];
1051: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1052: const PetscInt q = pind ? pind[p] : p;
1054: sf->atlasOff[q] = offset;
1055: offset += sf->atlasDof[q];
1056: }
1057: }
1058: /* Disable point offsets since these are unused */
1059: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1060: s->atlasOff[p] = -1;
1061: s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1062: }
1063: }
1064: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1065: /* Setup BC sections */
1066: PetscSectionSetUpBC(s);
1067: for (f = 0; f < s->numFields; ++f) {PetscSectionSetUpBC(s->field[f]);}
1068: return(0);
1069: }
1071: /*@
1072: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the chart
1074: Not collective
1076: Input Parameters:
1077: . s - the PetscSection
1079: Output Parameter:
1080: . maxDof - the maximum dof
1082: Level: intermediate
1084: .seealso: PetscSectionGetDof(), PetscSectionSetDof(), PetscSectionCreate()
1085: @*/
1086: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1087: {
1091: *maxDof = s->maxDof;
1092: return(0);
1093: }
1095: /*@
1096: PetscSectionGetStorageSize - Return the size of an array or local Vec capable of holding all the degrees of freedom.
1098: Not collective
1100: Input Parameter:
1101: . s - the PetscSection
1103: Output Parameter:
1104: . size - the size of an array which can hold all the dofs
1106: Level: intermediate
1108: .seealso: PetscSectionGetOffset(), PetscSectionGetConstrainedStorageSize(), PetscSectionCreate()
1109: @*/
1110: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1111: {
1112: PetscInt p, n = 0;
1117: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1118: *size = n;
1119: return(0);
1120: }
1122: /*@
1123: PetscSectionGetConstrainedStorageSize - Return the size of an array or local Vec capable of holding all unconstrained degrees of freedom.
1125: Not collective
1127: Input Parameters:
1128: + s - the PetscSection
1129: - point - the point
1131: Output Parameter:
1132: . size - the size of an array which can hold all unconstrained dofs
1134: Level: intermediate
1136: .seealso: PetscSectionGetStorageSize(), PetscSectionGetOffset(), PetscSectionCreate()
1137: @*/
1138: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1139: {
1140: PetscInt p, n = 0;
1145: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1146: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1147: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1148: }
1149: *size = n;
1150: return(0);
1151: }
1153: /*@
1154: PetscSectionCreateGlobalSection - Create a section describing the global field layout using
1155: the local section and an SF describing the section point overlap.
1157: Input Parameters:
1158: + s - The PetscSection for the local field layout
1159: . sf - The SF describing parallel layout of the section points (leaves are unowned local points)
1160: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1161: - localOffsets - If PETSC_TRUE, use local rather than global offsets for the points
1163: Output Parameter:
1164: . gsection - The PetscSection for the global field layout
1166: Note: This gives negative sizes and offsets to points not owned by this process
1168: Level: intermediate
1170: .seealso: PetscSectionCreate()
1171: @*/
1172: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1173: {
1174: PetscSection gs;
1175: const PetscInt *pind = NULL;
1176: PetscInt *recv = NULL, *neg = NULL;
1177: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1178: PetscErrorCode ierr;
1186: PetscSectionCreate(PetscObjectComm((PetscObject) s), &gs);
1187: PetscSectionGetChart(s, &pStart, &pEnd);
1188: PetscSectionSetChart(gs, pStart, pEnd);
1189: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1190: nlocal = nroots; /* The local/leaf space matches global/root space */
1191: /* Must allocate for all points visible to SF, which may be more than this section */
1192: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1193: PetscSFGetLeafRange(sf, NULL, &maxleaf);
1194: if (nroots < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %D < pEnd %D", nroots, pEnd);
1195: if (maxleaf >= nroots) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %D >= nroots %D", maxleaf, nroots);
1196: PetscMalloc2(nroots,&neg,nlocal,&recv);
1197: PetscArrayzero(neg,nroots);
1198: }
1199: /* Mark all local points with negative dof */
1200: for (p = pStart; p < pEnd; ++p) {
1201: PetscSectionGetDof(s, p, &dof);
1202: PetscSectionSetDof(gs, p, dof);
1203: PetscSectionGetConstraintDof(s, p, &cdof);
1204: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(gs, p, cdof);}
1205: if (neg) neg[p] = -(dof+1);
1206: }
1207: PetscSectionSetUpBC(gs);
1208: 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]);}
1209: if (nroots >= 0) {
1210: PetscArrayzero(recv,nlocal);
1211: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1212: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1213: for (p = pStart; p < pEnd; ++p) {
1214: if (recv[p] < 0) {
1215: gs->atlasDof[p-pStart] = recv[p];
1216: PetscSectionGetDof(s, p, &dof);
1217: 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);
1218: }
1219: }
1220: }
1221: /* Calculate new sizes, get process offset, and calculate point offsets */
1222: if (s->perm) {ISGetIndices(s->perm, &pind);}
1223: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1224: const PetscInt q = pind ? pind[p] : p;
1226: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1227: gs->atlasOff[q] = off;
1228: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q]-cdof : 0;
1229: }
1230: if (!localOffsets) {
1231: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1232: globalOff -= off;
1233: }
1234: for (p = pStart, off = 0; p < pEnd; ++p) {
1235: gs->atlasOff[p-pStart] += globalOff;
1236: if (neg) neg[p] = -(gs->atlasOff[p-pStart]+1);
1237: }
1238: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1239: /* Put in negative offsets for ghost points */
1240: if (nroots >= 0) {
1241: PetscArrayzero(recv,nlocal);
1242: PetscSFBcastBegin(sf, MPIU_INT, neg, recv);
1243: PetscSFBcastEnd(sf, MPIU_INT, neg, recv);
1244: for (p = pStart; p < pEnd; ++p) {
1245: if (recv[p] < 0) gs->atlasOff[p-pStart] = recv[p];
1246: }
1247: }
1248: PetscFree2(neg,recv);
1249: PetscSectionViewFromOptions(gs, NULL, "-global_section_view");
1250: *gsection = gs;
1251: return(0);
1252: }
1254: /*@
1255: PetscSectionCreateGlobalSectionCensored - Create a section describing the global field layout using
1256: the local section and an SF describing the section point overlap.
1258: Input Parameters:
1259: + s - The PetscSection for the local field layout
1260: . sf - The SF describing parallel layout of the section points
1261: . includeConstraints - By default this is PETSC_FALSE, meaning that the global field vector will not possess constrained dofs
1262: . numExcludes - The number of exclusion ranges
1263: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are numExcludes pairs
1265: Output Parameter:
1266: . gsection - The PetscSection for the global field layout
1268: Note: This gives negative sizes and offsets to points not owned by this process
1270: Level: advanced
1272: .seealso: PetscSectionCreate()
1273: @*/
1274: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1275: {
1276: const PetscInt *pind = NULL;
1277: PetscInt *neg = NULL, *tmpOff = NULL;
1278: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1279: PetscErrorCode ierr;
1285: PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);
1286: PetscSectionGetChart(s, &pStart, &pEnd);
1287: PetscSectionSetChart(*gsection, pStart, pEnd);
1288: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1289: if (nroots >= 0) {
1290: if (nroots < pEnd-pStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %D < %D section size", nroots, pEnd-pStart);
1291: PetscCalloc1(nroots, &neg);
1292: if (nroots > pEnd-pStart) {
1293: PetscCalloc1(nroots, &tmpOff);
1294: } else {
1295: tmpOff = &(*gsection)->atlasDof[-pStart];
1296: }
1297: }
1298: /* Mark ghost points with negative dof */
1299: for (p = pStart; p < pEnd; ++p) {
1300: for (e = 0; e < numExcludes; ++e) {
1301: if ((p >= excludes[e*2+0]) && (p < excludes[e*2+1])) {
1302: PetscSectionSetDof(*gsection, p, 0);
1303: break;
1304: }
1305: }
1306: if (e < numExcludes) continue;
1307: PetscSectionGetDof(s, p, &dof);
1308: PetscSectionSetDof(*gsection, p, dof);
1309: PetscSectionGetConstraintDof(s, p, &cdof);
1310: if (!includeConstraints && cdof > 0) {PetscSectionSetConstraintDof(*gsection, p, cdof);}
1311: if (neg) neg[p] = -(dof+1);
1312: }
1313: PetscSectionSetUpBC(*gsection);
1314: if (nroots >= 0) {
1315: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1316: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1317: if (nroots > pEnd - pStart) {
1318: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasDof[p-pStart] = tmpOff[p];}
1319: }
1320: }
1321: /* Calculate new sizes, get proccess offset, and calculate point offsets */
1322: if (s->perm) {ISGetIndices(s->perm, &pind);}
1323: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1324: const PetscInt q = pind ? pind[p] : p;
1326: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1327: (*gsection)->atlasOff[q] = off;
1328: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q]-cdof : 0;
1329: }
1330: MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));
1331: globalOff -= off;
1332: for (p = 0, off = 0; p < pEnd-pStart; ++p) {
1333: (*gsection)->atlasOff[p] += globalOff;
1334: if (neg) neg[p+pStart] = -((*gsection)->atlasOff[p]+1);
1335: }
1336: if (s->perm) {ISRestoreIndices(s->perm, &pind);}
1337: /* Put in negative offsets for ghost points */
1338: if (nroots >= 0) {
1339: if (nroots == pEnd-pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1340: PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff);
1341: PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff);
1342: if (nroots > pEnd - pStart) {
1343: for (p = pStart; p < pEnd; ++p) {if (tmpOff[p] < 0) (*gsection)->atlasOff[p-pStart] = tmpOff[p];}
1344: }
1345: }
1346: if (nroots >= 0 && nroots > pEnd-pStart) {PetscFree(tmpOff);}
1347: PetscFree(neg);
1348: return(0);
1349: }
1351: /*@
1352: PetscSectionGetPointLayout - Get the PetscLayout associated with the section points
1354: Collective on comm
1356: Input Parameters:
1357: + comm - The MPI_Comm
1358: - s - The PetscSection
1360: Output Parameter:
1361: . layout - The point layout for the section
1363: Note: This is usually caleld for the default global section.
1365: Level: advanced
1367: .seealso: PetscSectionGetValueLayout(), PetscSectionCreate()
1368: @*/
1369: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1370: {
1371: PetscInt pStart, pEnd, p, localSize = 0;
1375: PetscSectionGetChart(s, &pStart, &pEnd);
1376: for (p = pStart; p < pEnd; ++p) {
1377: PetscInt dof;
1379: PetscSectionGetDof(s, p, &dof);
1380: if (dof > 0) ++localSize;
1381: }
1382: PetscLayoutCreate(comm, layout);
1383: PetscLayoutSetLocalSize(*layout, localSize);
1384: PetscLayoutSetBlockSize(*layout, 1);
1385: PetscLayoutSetUp(*layout);
1386: return(0);
1387: }
1389: /*@
1390: PetscSectionGetValueLayout - Get the PetscLayout associated with the section dofs.
1392: Collective on comm
1394: Input Parameters:
1395: + comm - The MPI_Comm
1396: - s - The PetscSection
1398: Output Parameter:
1399: . layout - The dof layout for the section
1401: Note: This is usually called for the default global section.
1403: Level: advanced
1405: .seealso: PetscSectionGetPointLayout(), PetscSectionCreate()
1406: @*/
1407: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1408: {
1409: PetscInt pStart, pEnd, p, localSize = 0;
1415: PetscSectionGetChart(s, &pStart, &pEnd);
1416: for (p = pStart; p < pEnd; ++p) {
1417: PetscInt dof,cdof;
1419: PetscSectionGetDof(s, p, &dof);
1420: PetscSectionGetConstraintDof(s, p, &cdof);
1421: if (dof-cdof > 0) localSize += dof-cdof;
1422: }
1423: PetscLayoutCreate(comm, layout);
1424: PetscLayoutSetLocalSize(*layout, localSize);
1425: PetscLayoutSetBlockSize(*layout, 1);
1426: PetscLayoutSetUp(*layout);
1427: return(0);
1428: }
1430: /*@
1431: PetscSectionGetOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1433: Not collective
1435: Input Parameters:
1436: + s - the PetscSection
1437: - point - the point
1439: Output Parameter:
1440: . offset - the offset
1442: Level: intermediate
1444: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate()
1445: @*/
1446: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1447: {
1451: #if defined(PETSC_USE_DEBUG)
1452: 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);
1453: #endif
1454: *offset = s->atlasOff[point - s->pStart];
1455: return(0);
1456: }
1458: /*@
1459: PetscSectionSetOffset - Set the offset into an array or local Vec for the dof associated with the given point.
1461: Not collective
1463: Input Parameters:
1464: + s - the PetscSection
1465: . point - the point
1466: - offset - the offset
1468: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1470: Level: intermediate
1472: .seealso: PetscSectionGetFieldOffset(), PetscSectionCreate(), PetscSectionSetUp()
1473: @*/
1474: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1475: {
1478: 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);
1479: s->atlasOff[point - s->pStart] = offset;
1480: return(0);
1481: }
1483: /*@
1484: PetscSectionGetFieldOffset - Return the offset into an array or local Vec for the dof associated with the given point.
1486: Not collective
1488: Input Parameters:
1489: + s - the PetscSection
1490: . point - the point
1491: - field - the field
1493: Output Parameter:
1494: . offset - the offset
1496: Level: intermediate
1498: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1499: @*/
1500: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1501: {
1507: 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);
1508: PetscSectionGetOffset(s->field[field], point, offset);
1509: return(0);
1510: }
1512: /*@
1513: PetscSectionSetFieldOffset - Set 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
1520: . field - the field
1521: - offset - the offset
1523: Note: The user usually does not call this function, but uses PetscSectionSetUp()
1525: Level: intermediate
1527: .seealso: PetscSectionGetOffset(), PetscSectionCreate(), PetscSectionSetUp()
1528: @*/
1529: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1530: {
1535: 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);
1536: PetscSectionSetOffset(s->field[field], point, offset);
1537: return(0);
1538: }
1540: /*@
1541: PetscSectionGetFieldPointOffset - Return the offset on the given point for the dof associated with the given point.
1543: Not collective
1545: Input Parameters:
1546: + s - the PetscSection
1547: . point - the point
1548: - field - the field
1550: Output Parameter:
1551: . offset - the offset
1553: Note: This gives the offset on a point of the field, ignoring constraints, meaning starting at the first dof for
1554: this point, what number is the first dof with this field.
1556: Level: advanced
1558: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1559: @*/
1560: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1561: {
1562: PetscInt off, foff;
1568: 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);
1569: PetscSectionGetOffset(s, point, &off);
1570: PetscSectionGetOffset(s->field[field], point, &foff);
1571: *offset = foff - off;
1572: return(0);
1573: }
1575: /*@
1576: PetscSectionGetOffsetRange - Return the full range of offsets [start, end)
1578: Not collective
1580: Input Parameter:
1581: . s - the PetscSection
1583: Output Parameters:
1584: + start - the minimum offset
1585: - end - one more than the maximum offset
1587: Level: intermediate
1589: .seealso: PetscSectionGetOffset(), PetscSectionCreate()
1590: @*/
1591: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1592: {
1593: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1598: if (s->atlasOff) {os = s->atlasOff[0]; oe = s->atlasOff[0];}
1599: PetscSectionGetChart(s, &pStart, &pEnd);
1600: for (p = 0; p < pEnd-pStart; ++p) {
1601: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1603: if (off >= 0) {
1604: os = PetscMin(os, off);
1605: oe = PetscMax(oe, off+dof);
1606: }
1607: }
1608: if (start) *start = os;
1609: if (end) *end = oe;
1610: return(0);
1611: }
1613: /*@
1614: PetscSectionCreateSubsection - Create a new, smaller section composed of only the selected fields
1616: Collective on s
1618: Input Parameter:
1619: + s - the PetscSection
1620: . len - the number of subfields
1621: - fields - the subfield numbers
1623: Output Parameter:
1624: . subs - the subsection
1626: Note: The section offsets now refer to a new, smaller vector.
1628: Level: advanced
1630: .seealso: PetscSectionCreateSupersection(), PetscSectionCreate()
1631: @*/
1632: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
1633: {
1634: PetscInt nF, f, pStart, pEnd, p, maxCdof = 0;
1638: if (!len) return(0);
1642: PetscSectionGetNumFields(s, &nF);
1643: if (len > nF) SETERRQ2(PetscObjectComm((PetscObject) s), PETSC_ERR_ARG_WRONG, "Number of requested fields %D greater than number of fields %D", len, nF);
1644: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1645: PetscSectionSetNumFields(*subs, len);
1646: for (f = 0; f < len; ++f) {
1647: const char *name = NULL;
1648: PetscInt numComp = 0;
1650: PetscSectionGetFieldName(s, fields[f], &name);
1651: PetscSectionSetFieldName(*subs, f, name);
1652: PetscSectionGetFieldComponents(s, fields[f], &numComp);
1653: PetscSectionSetFieldComponents(*subs, f, numComp);
1654: }
1655: PetscSectionGetChart(s, &pStart, &pEnd);
1656: PetscSectionSetChart(*subs, pStart, pEnd);
1657: for (p = pStart; p < pEnd; ++p) {
1658: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
1660: for (f = 0; f < len; ++f) {
1661: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1662: PetscSectionSetFieldDof(*subs, p, f, fdof);
1663: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1664: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof);}
1665: dof += fdof;
1666: cdof += cfdof;
1667: }
1668: PetscSectionSetDof(*subs, p, dof);
1669: if (cdof) {PetscSectionSetConstraintDof(*subs, p, cdof);}
1670: maxCdof = PetscMax(cdof, maxCdof);
1671: }
1672: PetscSectionSetUp(*subs);
1673: if (maxCdof) {
1674: PetscInt *indices;
1676: PetscMalloc1(maxCdof, &indices);
1677: for (p = pStart; p < pEnd; ++p) {
1678: PetscInt cdof;
1680: PetscSectionGetConstraintDof(*subs, p, &cdof);
1681: if (cdof) {
1682: const PetscInt *oldIndices = NULL;
1683: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
1685: for (f = 0; f < len; ++f) {
1686: PetscInt oldFoff = 0, oldf;
1688: PetscSectionGetFieldDof(s, p, fields[f], &fdof);
1689: PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof);
1690: PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices);
1691: /* This can be sped up if we assume sorted fields */
1692: for (oldf = 0; oldf < fields[f]; ++oldf) {
1693: PetscInt oldfdof = 0;
1694: PetscSectionGetFieldDof(s, p, oldf, &oldfdof);
1695: oldFoff += oldfdof;
1696: }
1697: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + (fOff - oldFoff);
1698: PetscSectionSetFieldConstraintIndices(*subs, p, f, &indices[numConst]);
1699: numConst += cfdof;
1700: fOff += fdof;
1701: }
1702: PetscSectionSetConstraintIndices(*subs, p, indices);
1703: }
1704: }
1705: PetscFree(indices);
1706: }
1707: return(0);
1708: }
1710: /*@
1711: PetscSectionCreateSupersection - Create a new, larger section composed of the input sections
1713: Collective on s
1715: Input Parameters:
1716: + s - the input sections
1717: - len - the number of input sections
1719: Output Parameter:
1720: . supers - the supersection
1722: Note: The section offsets now refer to a new, larger vector.
1724: Level: advanced
1726: .seealso: PetscSectionCreateSubsection(), PetscSectionCreate()
1727: @*/
1728: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
1729: {
1730: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
1734: if (!len) return(0);
1735: for (i = 0; i < len; ++i) {
1736: PetscInt nf, pStarti, pEndi;
1738: PetscSectionGetNumFields(s[i], &nf);
1739: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1740: pStart = PetscMin(pStart, pStarti);
1741: pEnd = PetscMax(pEnd, pEndi);
1742: Nf += nf;
1743: }
1744: PetscSectionCreate(PetscObjectComm((PetscObject) s[0]), supers);
1745: PetscSectionSetNumFields(*supers, Nf);
1746: for (i = 0, f = 0; i < len; ++i) {
1747: PetscInt nf, fi;
1749: PetscSectionGetNumFields(s[i], &nf);
1750: for (fi = 0; fi < nf; ++fi, ++f) {
1751: const char *name = NULL;
1752: PetscInt numComp = 0;
1754: PetscSectionGetFieldName(s[i], fi, &name);
1755: PetscSectionSetFieldName(*supers, f, name);
1756: PetscSectionGetFieldComponents(s[i], fi, &numComp);
1757: PetscSectionSetFieldComponents(*supers, f, numComp);
1758: }
1759: }
1760: PetscSectionSetChart(*supers, pStart, pEnd);
1761: for (p = pStart; p < pEnd; ++p) {
1762: PetscInt dof = 0, cdof = 0;
1764: for (i = 0, f = 0; i < len; ++i) {
1765: PetscInt nf, fi, pStarti, pEndi;
1766: PetscInt fdof = 0, cfdof = 0;
1768: PetscSectionGetNumFields(s[i], &nf);
1769: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1770: if ((p < pStarti) || (p >= pEndi)) continue;
1771: for (fi = 0; fi < nf; ++fi, ++f) {
1772: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1773: PetscSectionAddFieldDof(*supers, p, f, fdof);
1774: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1775: if (cfdof) {PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof);}
1776: dof += fdof;
1777: cdof += cfdof;
1778: }
1779: }
1780: PetscSectionSetDof(*supers, p, dof);
1781: if (cdof) {PetscSectionSetConstraintDof(*supers, p, cdof);}
1782: maxCdof = PetscMax(cdof, maxCdof);
1783: }
1784: PetscSectionSetUp(*supers);
1785: if (maxCdof) {
1786: PetscInt *indices;
1788: PetscMalloc1(maxCdof, &indices);
1789: for (p = pStart; p < pEnd; ++p) {
1790: PetscInt cdof;
1792: PetscSectionGetConstraintDof(*supers, p, &cdof);
1793: if (cdof) {
1794: PetscInt dof, numConst = 0, fOff = 0;
1796: for (i = 0, f = 0; i < len; ++i) {
1797: const PetscInt *oldIndices = NULL;
1798: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
1800: PetscSectionGetNumFields(s[i], &nf);
1801: PetscSectionGetChart(s[i], &pStarti, &pEndi);
1802: if ((p < pStarti) || (p >= pEndi)) continue;
1803: for (fi = 0; fi < nf; ++fi, ++f) {
1804: PetscSectionGetFieldDof(s[i], p, fi, &fdof);
1805: PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof);
1806: PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices);
1807: for (fc = 0; fc < cfdof; ++fc) indices[numConst+fc] = oldIndices[fc] + fOff;
1808: PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]);
1809: numConst += cfdof;
1810: }
1811: PetscSectionGetDof(s[i], p, &dof);
1812: fOff += dof;
1813: }
1814: PetscSectionSetConstraintIndices(*supers, p, indices);
1815: }
1816: }
1817: PetscFree(indices);
1818: }
1819: return(0);
1820: }
1822: /*@
1823: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
1825: Collective on s
1827: Input Parameters:
1828: + s - the PetscSection
1829: - subpointMap - a sorted list of points in the original mesh which are in the submesh
1831: Output Parameter:
1832: . subs - the subsection
1834: Note: The section offsets now refer to a new, smaller vector.
1836: Level: advanced
1838: .seealso: PetscSectionCreateSubsection(), DMPlexGetSubpointMap(), PetscSectionCreate()
1839: @*/
1840: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
1841: {
1842: const PetscInt *points = NULL, *indices = NULL;
1843: PetscInt numFields, f, numSubpoints = 0, pStart, pEnd, p, subp;
1850: PetscSectionGetNumFields(s, &numFields);
1851: PetscSectionCreate(PetscObjectComm((PetscObject) s), subs);
1852: if (numFields) {PetscSectionSetNumFields(*subs, numFields);}
1853: for (f = 0; f < numFields; ++f) {
1854: const char *name = NULL;
1855: PetscInt numComp = 0;
1857: PetscSectionGetFieldName(s, f, &name);
1858: PetscSectionSetFieldName(*subs, f, name);
1859: PetscSectionGetFieldComponents(s, f, &numComp);
1860: PetscSectionSetFieldComponents(*subs, f, numComp);
1861: }
1862: /* For right now, we do not try to squeeze the subchart */
1863: if (subpointMap) {
1864: ISGetSize(subpointMap, &numSubpoints);
1865: ISGetIndices(subpointMap, &points);
1866: }
1867: PetscSectionGetChart(s, &pStart, &pEnd);
1868: PetscSectionSetChart(*subs, 0, numSubpoints);
1869: for (p = pStart; p < pEnd; ++p) {
1870: PetscInt dof, cdof, fdof = 0, cfdof = 0;
1872: PetscFindInt(p, numSubpoints, points, &subp);
1873: if (subp < 0) continue;
1874: for (f = 0; f < numFields; ++f) {
1875: PetscSectionGetFieldDof(s, p, f, &fdof);
1876: PetscSectionSetFieldDof(*subs, subp, f, fdof);
1877: PetscSectionGetFieldConstraintDof(s, p, f, &cfdof);
1878: if (cfdof) {PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof);}
1879: }
1880: PetscSectionGetDof(s, p, &dof);
1881: PetscSectionSetDof(*subs, subp, dof);
1882: PetscSectionGetConstraintDof(s, p, &cdof);
1883: if (cdof) {PetscSectionSetConstraintDof(*subs, subp, cdof);}
1884: }
1885: PetscSectionSetUp(*subs);
1886: /* Change offsets to original offsets */
1887: for (p = pStart; p < pEnd; ++p) {
1888: PetscInt off, foff = 0;
1890: PetscFindInt(p, numSubpoints, points, &subp);
1891: if (subp < 0) continue;
1892: for (f = 0; f < numFields; ++f) {
1893: PetscSectionGetFieldOffset(s, p, f, &foff);
1894: PetscSectionSetFieldOffset(*subs, subp, f, foff);
1895: }
1896: PetscSectionGetOffset(s, p, &off);
1897: PetscSectionSetOffset(*subs, subp, off);
1898: }
1899: /* Copy constraint indices */
1900: for (subp = 0; subp < numSubpoints; ++subp) {
1901: PetscInt cdof;
1903: PetscSectionGetConstraintDof(*subs, subp, &cdof);
1904: if (cdof) {
1905: for (f = 0; f < numFields; ++f) {
1906: PetscSectionGetFieldConstraintIndices(s, points[subp], f, &indices);
1907: PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices);
1908: }
1909: PetscSectionGetConstraintIndices(s, points[subp], &indices);
1910: PetscSectionSetConstraintIndices(*subs, subp, indices);
1911: }
1912: }
1913: if (subpointMap) {ISRestoreIndices(subpointMap, &points);}
1914: return(0);
1915: }
1917: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
1918: {
1919: PetscInt p;
1920: PetscMPIInt rank;
1924: MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
1925: PetscViewerASCIIPushSynchronized(viewer);
1926: PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank);
1927: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1928: if ((s->bc) && (s->bc->atlasDof[p] > 0)) {
1929: PetscInt b;
1931: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D constrained", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1932: for (b = 0; b < s->bc->atlasDof[p]; ++b) {
1933: PetscViewerASCIISynchronizedPrintf(viewer, " %D", s->bcIndices[s->bc->atlasOff[p]+b]);
1934: }
1935: PetscViewerASCIISynchronizedPrintf(viewer, "\n");
1936: } else {
1937: PetscViewerASCIISynchronizedPrintf(viewer, " (%4D) dim %2D offset %3D\n", p+s->pStart, s->atlasDof[p], s->atlasOff[p]);
1938: }
1939: }
1940: PetscViewerFlush(viewer);
1941: PetscViewerASCIIPopSynchronized(viewer);
1942: if (s->sym) {
1943: PetscViewerASCIIPushTab(viewer);
1944: PetscSectionSymView(s->sym,viewer);
1945: PetscViewerASCIIPopTab(viewer);
1946: }
1947: return(0);
1948: }
1950: /*@C
1951: PetscSectionView - Views a PetscSection
1953: Collective on PetscSection
1955: Input Parameters:
1956: + s - the PetscSection object to view
1957: - v - the viewer
1959: Level: beginner
1961: .seealso PetscSectionCreate(), PetscSectionDestroy()
1962: @*/
1963: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
1964: {
1965: PetscBool isascii;
1966: PetscInt f;
1971: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer);}
1973: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1974: if (isascii) {
1975: PetscObjectPrintClassNamePrefixType((PetscObject)s,viewer);
1976: if (s->numFields) {
1977: PetscViewerASCIIPrintf(viewer, "%D fields\n", s->numFields);
1978: for (f = 0; f < s->numFields; ++f) {
1979: PetscViewerASCIIPrintf(viewer, " field %D with %D components\n", f, s->numFieldComponents[f]);
1980: PetscSectionView_ASCII(s->field[f], viewer);
1981: }
1982: } else {
1983: PetscSectionView_ASCII(s, viewer);
1984: }
1985: }
1986: return(0);
1987: }
1989: /*@
1990: PetscSectionReset - Frees all section data.
1992: Not collective
1994: Input Parameters:
1995: . s - the PetscSection
1997: Level: beginner
1999: .seealso: PetscSection, PetscSectionCreate()
2000: @*/
2001: PetscErrorCode PetscSectionReset(PetscSection s)
2002: {
2003: PetscInt f;
2008: PetscFree(s->numFieldComponents);
2009: for (f = 0; f < s->numFields; ++f) {
2010: PetscSectionDestroy(&s->field[f]);
2011: PetscFree(s->fieldNames[f]);
2012: }
2013: PetscFree(s->fieldNames);
2014: PetscFree(s->field);
2015: PetscSectionDestroy(&s->bc);
2016: PetscFree(s->bcIndices);
2017: PetscFree2(s->atlasDof, s->atlasOff);
2018: PetscSectionDestroy(&s->clSection);
2019: ISDestroy(&s->clPoints);
2020: ISDestroy(&s->perm);
2021: PetscFree(s->clPerm);
2022: PetscFree(s->clInvPerm);
2023: PetscSectionSymDestroy(&s->sym);
2025: s->pStart = -1;
2026: s->pEnd = -1;
2027: s->maxDof = 0;
2028: s->setup = PETSC_FALSE;
2029: s->numFields = 0;
2030: s->clObj = NULL;
2031: return(0);
2032: }
2034: /*@
2035: PetscSectionDestroy - Frees a section object and frees its range if that exists.
2037: Not collective
2039: Input Parameters:
2040: . s - the PetscSection
2042: Level: beginner
2044: .seealso: PetscSection, PetscSectionCreate()
2045: @*/
2046: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2047: {
2051: if (!*s) return(0);
2053: if (--((PetscObject)(*s))->refct > 0) {
2054: *s = NULL;
2055: return(0);
2056: }
2057: PetscSectionReset(*s);
2058: PetscHeaderDestroy(s);
2059: return(0);
2060: }
2062: PetscErrorCode VecIntGetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2063: {
2064: const PetscInt p = point - s->pStart;
2068: *values = &baseArray[s->atlasOff[p]];
2069: return(0);
2070: }
2072: PetscErrorCode VecIntSetValuesSection(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2073: {
2074: PetscInt *array;
2075: const PetscInt p = point - s->pStart;
2076: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2077: PetscInt cDim = 0;
2082: PetscSectionGetConstraintDof(s, p, &cDim);
2083: array = &baseArray[s->atlasOff[p]];
2084: if (!cDim) {
2085: if (orientation >= 0) {
2086: const PetscInt dim = s->atlasDof[p];
2087: PetscInt i;
2089: if (mode == INSERT_VALUES) {
2090: for (i = 0; i < dim; ++i) array[i] = values[i];
2091: } else {
2092: for (i = 0; i < dim; ++i) array[i] += values[i];
2093: }
2094: } else {
2095: PetscInt offset = 0;
2096: PetscInt j = -1, field, i;
2098: for (field = 0; field < s->numFields; ++field) {
2099: const PetscInt dim = s->field[field]->atlasDof[p];
2101: for (i = dim-1; i >= 0; --i) array[++j] = values[i+offset];
2102: offset += dim;
2103: }
2104: }
2105: } else {
2106: if (orientation >= 0) {
2107: const PetscInt dim = s->atlasDof[p];
2108: PetscInt cInd = 0, i;
2109: const PetscInt *cDof;
2111: PetscSectionGetConstraintIndices(s, point, &cDof);
2112: if (mode == INSERT_VALUES) {
2113: for (i = 0; i < dim; ++i) {
2114: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2115: array[i] = values[i];
2116: }
2117: } else {
2118: for (i = 0; i < dim; ++i) {
2119: if ((cInd < cDim) && (i == cDof[cInd])) {++cInd; continue;}
2120: array[i] += values[i];
2121: }
2122: }
2123: } else {
2124: const PetscInt *cDof;
2125: PetscInt offset = 0;
2126: PetscInt cOffset = 0;
2127: PetscInt j = 0, field;
2129: PetscSectionGetConstraintIndices(s, point, &cDof);
2130: for (field = 0; field < s->numFields; ++field) {
2131: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2132: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2133: const PetscInt sDim = dim - tDim;
2134: PetscInt cInd = 0, i,k;
2136: for (i = 0, k = dim+offset-1; i < dim; ++i, ++j, --k) {
2137: if ((cInd < sDim) && (j == cDof[cInd+cOffset])) {++cInd; continue;}
2138: array[j] = values[k];
2139: }
2140: offset += dim;
2141: cOffset += dim - tDim;
2142: }
2143: }
2144: }
2145: return(0);
2146: }
2148: /*@C
2149: PetscSectionHasConstraints - Determine whether a section has constrained dofs
2151: Not collective
2153: Input Parameter:
2154: . s - The PetscSection
2156: Output Parameter:
2157: . hasConstraints - flag indicating that the section has constrained dofs
2159: Level: intermediate
2161: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2162: @*/
2163: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2164: {
2168: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2169: return(0);
2170: }
2172: /*@C
2173: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained
2175: Not collective
2177: Input Parameters:
2178: + s - The PetscSection
2179: - point - The point
2181: Output Parameter:
2182: . indices - The constrained dofs
2184: Note: In Fortran, you call PetscSectionGetConstraintIndicesF90() and PetscSectionRestoreConstraintIndicesF90()
2186: Level: intermediate
2188: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2189: @*/
2190: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2191: {
2196: if (s->bc) {
2197: VecIntGetValuesSection(s->bcIndices, s->bc, point, indices);
2198: } else *indices = NULL;
2199: return(0);
2200: }
2202: /*@C
2203: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2205: Not collective
2207: Input Parameters:
2208: + s - The PetscSection
2209: . point - The point
2210: - indices - The constrained dofs
2212: Note: The Fortran is PetscSectionSetConstraintIndicesF90()
2214: Level: intermediate
2216: .seealso: PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2217: @*/
2218: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2219: {
2224: if (s->bc) {
2225: VecIntSetValuesSection(s->bcIndices, s->bc, point, indices, INSERT_VALUES);
2226: }
2227: return(0);
2228: }
2230: /*@C
2231: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2233: Not collective
2235: Input Parameters:
2236: + s - The PetscSection
2237: . field - The field number
2238: - point - The point
2240: Output Parameter:
2241: . indices - The constrained dofs
2243: Note: In Fortran, you call PetscSectionGetFieldConstraintIndicesF90() and PetscSectionRestoreFieldConstraintIndicesF90()
2245: Level: intermediate
2247: .seealso: PetscSectionSetFieldConstraintIndices(), PetscSectionGetConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2248: @*/
2249: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2250: {
2255: 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);
2256: PetscSectionGetConstraintIndices(s->field[field], point, indices);
2257: return(0);
2258: }
2260: /*@C
2261: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2263: Not collective
2265: Input Parameters:
2266: + s - The PetscSection
2267: . point - The point
2268: . field - The field number
2269: - indices - The constrained dofs
2271: Note: The Fortran is PetscSectionSetFieldConstraintIndicesF90()
2273: Level: intermediate
2275: .seealso: PetscSectionSetConstraintIndices(), PetscSectionGetFieldConstraintIndices(), PetscSectionGetConstraintDof(), PetscSection
2276: @*/
2277: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2278: {
2283: 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);
2284: PetscSectionSetConstraintIndices(s->field[field], point, indices);
2285: return(0);
2286: }
2288: /*@
2289: PetscSectionPermute - Reorder the section according to the input point permutation
2291: Collective on PetscSection
2293: Input Parameter:
2294: + section - The PetscSection object
2295: - perm - The point permutation, old point p becomes new point perm[p]
2297: Output Parameter:
2298: . sectionNew - The permuted PetscSection
2300: Level: intermediate
2302: .seealso: MatPermute()
2303: @*/
2304: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2305: {
2306: PetscSection s = section, sNew;
2307: const PetscInt *perm;
2308: PetscInt numFields, f, numPoints, pStart, pEnd, p;
2309: PetscErrorCode ierr;
2315: PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);
2316: PetscSectionGetNumFields(s, &numFields);
2317: if (numFields) {PetscSectionSetNumFields(sNew, numFields);}
2318: for (f = 0; f < numFields; ++f) {
2319: const char *name;
2320: PetscInt numComp;
2322: PetscSectionGetFieldName(s, f, &name);
2323: PetscSectionSetFieldName(sNew, f, name);
2324: PetscSectionGetFieldComponents(s, f, &numComp);
2325: PetscSectionSetFieldComponents(sNew, f, numComp);
2326: }
2327: ISGetLocalSize(permutation, &numPoints);
2328: ISGetIndices(permutation, &perm);
2329: PetscSectionGetChart(s, &pStart, &pEnd);
2330: PetscSectionSetChart(sNew, pStart, pEnd);
2331: if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %D is less than largest Section point %D", numPoints, pEnd);
2332: for (p = pStart; p < pEnd; ++p) {
2333: PetscInt dof, cdof;
2335: PetscSectionGetDof(s, p, &dof);
2336: PetscSectionSetDof(sNew, perm[p], dof);
2337: PetscSectionGetConstraintDof(s, p, &cdof);
2338: if (cdof) {PetscSectionSetConstraintDof(sNew, perm[p], cdof);}
2339: for (f = 0; f < numFields; ++f) {
2340: PetscSectionGetFieldDof(s, p, f, &dof);
2341: PetscSectionSetFieldDof(sNew, perm[p], f, dof);
2342: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2343: if (cdof) {PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);}
2344: }
2345: }
2346: PetscSectionSetUp(sNew);
2347: for (p = pStart; p < pEnd; ++p) {
2348: const PetscInt *cind;
2349: PetscInt cdof;
2351: PetscSectionGetConstraintDof(s, p, &cdof);
2352: if (cdof) {
2353: PetscSectionGetConstraintIndices(s, p, &cind);
2354: PetscSectionSetConstraintIndices(sNew, perm[p], cind);
2355: }
2356: for (f = 0; f < numFields; ++f) {
2357: PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
2358: if (cdof) {
2359: PetscSectionGetFieldConstraintIndices(s, p, f, &cind);
2360: PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);
2361: }
2362: }
2363: }
2364: ISRestoreIndices(permutation, &perm);
2365: *sectionNew = sNew;
2366: return(0);
2367: }
2369: /* TODO: the next three functions should be moved to sf/utils */
2370: #include <petsc/private/sfimpl.h>
2372: /*@C
2373: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
2375: Collective on sf
2377: Input Parameters:
2378: + sf - The SF
2379: - rootSection - Section defined on root space
2381: Output Parameters:
2382: + remoteOffsets - root offsets in leaf storage, or NULL
2383: - leafSection - Section defined on the leaf space
2385: Level: advanced
2387: .seealso: PetscSFCreate()
2388: @*/
2389: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
2390: {
2391: PetscSF embedSF;
2392: const PetscInt *indices;
2393: IS selected;
2394: PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f;
2395: PetscBool *sub, hasc;
2399: PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
2400: PetscSectionGetNumFields(rootSection, &numFields);
2401: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
2402: PetscMalloc1(numFields+2, &sub);
2403: sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
2404: for (f = 0; f < numFields; ++f) {
2405: PetscSectionSym sym;
2406: const char *name = NULL;
2407: PetscInt numComp = 0;
2409: sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
2410: PetscSectionGetFieldComponents(rootSection, f, &numComp);
2411: PetscSectionGetFieldName(rootSection, f, &name);
2412: PetscSectionGetFieldSym(rootSection, f, &sym);
2413: PetscSectionSetFieldComponents(leafSection, f, numComp);
2414: PetscSectionSetFieldName(leafSection, f, name);
2415: PetscSectionSetFieldSym(leafSection, f, sym);
2416: }
2417: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2418: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
2419: rpEnd = PetscMin(rpEnd,nroots);
2420: rpEnd = PetscMax(rpStart,rpEnd);
2421: /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
2422: sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
2423: MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
2424: if (sub[0]) {
2425: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2426: ISGetIndices(selected, &indices);
2427: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2428: ISRestoreIndices(selected, &indices);
2429: ISDestroy(&selected);
2430: } else {
2431: PetscObjectReference((PetscObject)sf);
2432: embedSF = sf;
2433: }
2434: PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
2435: lpEnd++;
2437: PetscSectionSetChart(leafSection, lpStart, lpEnd);
2439: /* Constrained dof section */
2440: hasc = sub[1];
2441: for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);
2443: /* Could fuse these at the cost of copies and extra allocation */
2444: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2445: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart]);
2446: if (sub[1]) {
2447: PetscSectionCheckConstraints_Static(rootSection);
2448: PetscSectionCheckConstraints_Static(leafSection);
2449: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2450: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart]);
2451: }
2452: for (f = 0; f < numFields; ++f) {
2453: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2454: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart]);
2455: if (sub[2+f]) {
2456: PetscSectionCheckConstraints_Static(rootSection->field[f]);
2457: PetscSectionCheckConstraints_Static(leafSection->field[f]);
2458: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2459: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart]);
2460: }
2461: }
2462: if (remoteOffsets) {
2463: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
2464: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2465: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2466: }
2467: PetscSectionSetUp(leafSection);
2468: if (hasc) { /* need to communicate bcIndices */
2469: PetscSF bcSF;
2470: PetscInt *rOffBc;
2472: PetscMalloc1(lpEnd - lpStart, &rOffBc);
2473: if (sub[1]) {
2474: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2475: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2476: PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
2477: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2478: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices);
2479: PetscSFDestroy(&bcSF);
2480: }
2481: for (f = 0; f < numFields; ++f) {
2482: if (sub[2+f]) {
2483: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2484: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart]);
2485: PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
2486: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2487: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices);
2488: PetscSFDestroy(&bcSF);
2489: }
2490: }
2491: PetscFree(rOffBc);
2492: }
2493: PetscSFDestroy(&embedSF);
2494: PetscFree(sub);
2495: PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
2496: return(0);
2497: }
2499: /*@C
2500: PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
2502: Collective on sf
2504: Input Parameters:
2505: + sf - The SF
2506: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
2507: - leafSection - Data layout of local points for incoming data (this is layout for SF leaves)
2509: Output Parameter:
2510: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2512: Level: developer
2514: .seealso: PetscSFCreate()
2515: @*/
2516: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
2517: {
2518: PetscSF embedSF;
2519: const PetscInt *indices;
2520: IS selected;
2521: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
2522: PetscErrorCode ierr;
2525: *remoteOffsets = NULL;
2526: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
2527: if (numRoots < 0) return(0);
2528: PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
2529: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
2530: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2531: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
2532: ISGetIndices(selected, &indices);
2533: PetscSFCreateEmbeddedSF(sf, rpEnd - rpStart, indices, &embedSF);
2534: ISRestoreIndices(selected, &indices);
2535: ISDestroy(&selected);
2536: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
2537: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2538: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart]);
2539: PetscSFDestroy(&embedSF);
2540: PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
2541: return(0);
2542: }
2544: /*@C
2545: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
2547: Collective on sf
2549: Input Parameters:
2550: + sf - The SF
2551: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
2552: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
2553: - leafSection - Data layout of local points for incoming data (this is the distributed section)
2555: Output Parameters:
2556: - sectionSF - The new SF
2558: Note: Either rootSection or remoteOffsets can be specified
2560: Level: advanced
2562: .seealso: PetscSFCreate()
2563: @*/
2564: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2565: {
2566: MPI_Comm comm;
2567: const PetscInt *localPoints;
2568: const PetscSFNode *remotePoints;
2569: PetscInt lpStart, lpEnd;
2570: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
2571: PetscInt *localIndices;
2572: PetscSFNode *remoteIndices;
2573: PetscInt i, ind;
2574: PetscErrorCode ierr;
2582: PetscObjectGetComm((PetscObject)sf,&comm);
2583: PetscSFCreate(comm, sectionSF);
2584: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
2585: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
2586: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
2587: if (numRoots < 0) return(0);
2588: PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
2589: for (i = 0; i < numPoints; ++i) {
2590: PetscInt localPoint = localPoints ? localPoints[i] : i;
2591: PetscInt dof;
2593: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2594: PetscSectionGetDof(leafSection, localPoint, &dof);
2595: numIndices += dof;
2596: }
2597: }
2598: PetscMalloc1(numIndices, &localIndices);
2599: PetscMalloc1(numIndices, &remoteIndices);
2600: /* Create new index graph */
2601: for (i = 0, ind = 0; i < numPoints; ++i) {
2602: PetscInt localPoint = localPoints ? localPoints[i] : i;
2603: PetscInt rank = remotePoints[i].rank;
2605: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
2606: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
2607: PetscInt loff, dof, d;
2609: PetscSectionGetOffset(leafSection, localPoint, &loff);
2610: PetscSectionGetDof(leafSection, localPoint, &dof);
2611: for (d = 0; d < dof; ++d, ++ind) {
2612: localIndices[ind] = loff+d;
2613: remoteIndices[ind].rank = rank;
2614: remoteIndices[ind].index = remoteOffset+d;
2615: }
2616: }
2617: }
2618: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
2619: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
2620: PetscSFSetUp(*sectionSF);
2621: PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
2622: return(0);
2623: }
2625: /*@
2626: PetscSectionSetClosureIndex - Set a cache of points in the closure of each point in the section
2628: Collective on section
2630: Input Parameters:
2631: + section - The PetscSection
2632: . obj - A PetscObject which serves as the key for this index
2633: . clSection - Section giving the size of the closure of each point
2634: - clPoints - IS giving the points in each closure
2636: Note: We compress out closure points with no dofs in this section
2638: Level: advanced
2640: .seealso: PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex()
2641: @*/
2642: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
2643: {
2647: if (section->clObj != obj) {PetscFree(section->clPerm);PetscFree(section->clInvPerm);}
2648: section->clObj = obj;
2649: PetscSectionDestroy(§ion->clSection);
2650: ISDestroy(§ion->clPoints);
2651: section->clSection = clSection;
2652: section->clPoints = clPoints;
2653: return(0);
2654: }
2656: /*@
2657: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section
2659: Collective on section
2661: Input Parameters:
2662: + section - The PetscSection
2663: - obj - A PetscObject which serves as the key for this index
2665: Output Parameters:
2666: + clSection - Section giving the size of the closure of each point
2667: - clPoints - IS giving the points in each closure
2669: Note: We compress out closure points with no dofs in this section
2671: Level: advanced
2673: .seealso: PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2674: @*/
2675: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
2676: {
2678: if (section->clObj == obj) {
2679: if (clSection) *clSection = section->clSection;
2680: if (clPoints) *clPoints = section->clPoints;
2681: } else {
2682: if (clSection) *clSection = NULL;
2683: if (clPoints) *clPoints = NULL;
2684: }
2685: return(0);
2686: }
2688: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
2689: {
2690: PetscInt i;
2694: if (section->clObj != obj) {
2695: PetscSectionDestroy(§ion->clSection);
2696: ISDestroy(§ion->clPoints);
2697: }
2698: section->clObj = obj;
2699: PetscFree(section->clPerm);
2700: PetscFree(section->clInvPerm);
2701: section->clSize = clSize;
2702: if (mode == PETSC_COPY_VALUES) {
2703: PetscMalloc1(clSize, §ion->clPerm);
2704: PetscLogObjectMemory((PetscObject) obj, clSize*sizeof(PetscInt));
2705: PetscArraycpy(section->clPerm, clPerm, clSize);
2706: } else if (mode == PETSC_OWN_POINTER) {
2707: section->clPerm = clPerm;
2708: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
2709: PetscMalloc1(clSize, §ion->clInvPerm);
2710: for (i = 0; i < clSize; ++i) section->clInvPerm[section->clPerm[i]] = i;
2711: return(0);
2712: }
2714: /*@
2715: PetscSectionSetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2717: Not Collective
2719: Input Parameters:
2720: + section - The PetscSection
2721: . obj - A PetscObject which serves as the key for this index
2722: - perm - Permutation of the cell dof closure
2724: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2725: other points (like faces).
2727: Level: intermediate
2729: .seealso: PetscSectionGetClosurePermutation(), PetscSectionGetClosureIndex(), DMPlexCreateClosureIndex(), PetscCopyMode
2730: @*/
2731: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, IS perm)
2732: {
2733: const PetscInt *clPerm = NULL;
2734: PetscInt clSize = 0;
2735: PetscErrorCode ierr;
2738: if (perm) {
2739: ISGetLocalSize(perm, &clSize);
2740: ISGetIndices(perm, &clPerm);
2741: }
2742: PetscSectionSetClosurePermutation_Internal(section, obj, clSize, PETSC_COPY_VALUES, (PetscInt *) clPerm);
2743: if (perm) {ISRestoreIndices(perm, &clPerm);}
2744: return(0);
2745: }
2747: PetscErrorCode PetscSectionGetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2748: {
2750: if (section->clObj == obj) {
2751: if (size) *size = section->clSize;
2752: if (perm) *perm = section->clPerm;
2753: } else {
2754: if (size) *size = 0;
2755: if (perm) *perm = NULL;
2756: }
2757: return(0);
2758: }
2760: /*@
2761: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
2763: Not collective
2765: Input Parameters:
2766: + section - The PetscSection
2767: - obj - A PetscObject which serves as the key for this index
2769: Output Parameter:
2770: . perm - The dof closure permutation
2772: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2773: other points (like faces).
2775: The user must destroy the IS that is returned.
2777: Level: intermediate
2779: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureInversePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2780: @*/
2781: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, IS *perm)
2782: {
2783: const PetscInt *clPerm;
2784: PetscInt clSize;
2785: PetscErrorCode ierr;
2788: PetscSectionGetClosurePermutation_Internal(section, obj, &clSize, &clPerm);
2789: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2790: return(0);
2791: }
2793: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt *size, const PetscInt *perm[])
2794: {
2796: if (section->clObj == obj) {
2797: if (size) *size = section->clSize;
2798: if (perm) *perm = section->clInvPerm;
2799: } else {
2800: if (size) *size = 0;
2801: if (perm) *perm = NULL;
2802: }
2803: return(0);
2804: }
2806: /*@
2807: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
2809: Not collective
2811: Input Parameters:
2812: + section - The PetscSection
2813: - obj - A PetscObject which serves as the key for this index
2815: Output Parameters:
2816: + size - The dof closure size
2817: - perm - The dof closure permutation
2819: Note: This strategy only works when all cells have the same size dof closure, and no closures are retrieved for
2820: other points (like faces).
2822: The user must destroy the IS that is returned.
2824: Level: intermediate
2826: .seealso: PetscSectionSetClosurePermutation(), PetscSectionGetClosureIndex(), PetscSectionSetClosureIndex(), DMPlexCreateClosureIndex()
2827: @*/
2828: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, IS *perm)
2829: {
2830: const PetscInt *clPerm;
2831: PetscInt clSize;
2832: PetscErrorCode ierr;
2835: PetscSectionGetClosureInversePermutation_Internal(section, obj, &clSize, &clPerm);
2836: ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm);
2837: return(0);
2838: }
2840: /*@
2841: PetscSectionGetField - Get the subsection associated with a single field
2843: Input Parameters:
2844: + s - The PetscSection
2845: - field - The field number
2847: Output Parameter:
2848: . subs - The subsection for the given field
2850: Level: intermediate
2852: .seealso: PetscSectionSetNumFields()
2853: @*/
2854: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
2855: {
2859: 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);
2860: *subs = s->field[field];
2861: return(0);
2862: }
2864: PetscClassId PETSC_SECTION_SYM_CLASSID;
2865: PetscFunctionList PetscSectionSymList = NULL;
2867: /*@
2868: PetscSectionSymCreate - Creates an empty PetscSectionSym object.
2870: Collective
2872: Input Parameter:
2873: . comm - the MPI communicator
2875: Output Parameter:
2876: . sym - pointer to the new set of symmetries
2878: Level: developer
2880: .seealso: PetscSectionSym, PetscSectionSymDestroy()
2881: @*/
2882: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
2883: {
2888: ISInitializePackage();
2889: PetscHeaderCreate(*sym,PETSC_SECTION_SYM_CLASSID,"PetscSectionSym","Section Symmetry","IS",comm,PetscSectionSymDestroy,PetscSectionSymView);
2890: return(0);
2891: }
2893: /*@C
2894: PetscSectionSymSetType - Builds a PetscSection symmetry, for a particular implementation.
2896: Collective on PetscSectionSym
2898: Input Parameters:
2899: + sym - The section symmetry object
2900: - method - The name of the section symmetry type
2902: Level: developer
2904: .seealso: PetscSectionSymGetType(), PetscSectionSymCreate()
2905: @*/
2906: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
2907: {
2908: PetscErrorCode (*r)(PetscSectionSym);
2909: PetscBool match;
2914: PetscObjectTypeCompare((PetscObject) sym, method, &match);
2915: if (match) return(0);
2917: PetscFunctionListFind(PetscSectionSymList,method,&r);
2918: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
2919: if (sym->ops->destroy) {
2920: (*sym->ops->destroy)(sym);
2921: sym->ops->destroy = NULL;
2922: }
2923: (*r)(sym);
2924: PetscObjectChangeTypeName((PetscObject)sym,method);
2925: return(0);
2926: }
2929: /*@C
2930: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the PetscSectionSym.
2932: Not Collective
2934: Input Parameter:
2935: . sym - The section symmetry
2937: Output Parameter:
2938: . type - The index set type name
2940: Level: developer
2942: .seealso: PetscSectionSymSetType(), PetscSectionSymCreate()
2943: @*/
2944: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
2945: {
2949: *type = ((PetscObject)sym)->type_name;
2950: return(0);
2951: }
2953: /*@C
2954: PetscSectionSymRegister - Adds a new section symmetry implementation
2956: Not Collective
2958: Input Parameters:
2959: + name - The name of a new user-defined creation routine
2960: - create_func - The creation routine itself
2962: Notes:
2963: PetscSectionSymRegister() may be called multiple times to add several user-defined vectors
2965: Level: developer
2967: .seealso: PetscSectionSymCreate(), PetscSectionSymSetType()
2968: @*/
2969: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
2970: {
2974: ISInitializePackage();
2975: PetscFunctionListAdd(&PetscSectionSymList,sname,function);
2976: return(0);
2977: }
2979: /*@
2980: PetscSectionSymDestroy - Destroys a section symmetry.
2982: Collective on PetscSectionSym
2984: Input Parameters:
2985: . sym - the section symmetry
2987: Level: developer
2989: .seealso: PetscSectionSymCreate(), PetscSectionSymDestroy()
2990: @*/
2991: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
2992: {
2993: SymWorkLink link,next;
2997: if (!*sym) return(0);
2999: if (--((PetscObject)(*sym))->refct > 0) {*sym = 0; return(0);}
3000: if ((*sym)->ops->destroy) {
3001: (*(*sym)->ops->destroy)(*sym);
3002: }
3003: if ((*sym)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
3004: for (link=(*sym)->workin; link; link=next) {
3005: next = link->next;
3006: PetscFree2(*(PetscInt***)&link->perms,*(PetscScalar***)&link->rots);
3007: PetscFree(link);
3008: }
3009: (*sym)->workin = NULL;
3010: PetscHeaderDestroy(sym);
3011: return(0);
3012: }
3014: /*@C
3015: PetscSectionSymView - Displays a section symmetry
3017: Collective on PetscSectionSym
3019: Input Parameters:
3020: + sym - the index set
3021: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
3023: Level: developer
3025: .seealso: PetscViewerASCIIOpen()
3026: @*/
3027: PetscErrorCode PetscSectionSymView(PetscSectionSym sym,PetscViewer viewer)
3028: {
3033: if (!viewer) {
3034: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym),&viewer);
3035: }
3038: PetscObjectPrintClassNamePrefixType((PetscObject)sym,viewer);
3039: if (sym->ops->view) {
3040: (*sym->ops->view)(sym,viewer);
3041: }
3042: return(0);
3043: }
3045: /*@
3046: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3048: Collective on PetscSection
3050: Input Parameters:
3051: + section - the section describing data layout
3052: - sym - the symmetry describing the affect of orientation on the access of the data
3054: Level: developer
3056: .seealso: PetscSectionGetSym(), PetscSectionSymCreate()
3057: @*/
3058: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3059: {
3064: PetscSectionSymDestroy(&(section->sym));
3065: if (sym) {
3068: PetscObjectReference((PetscObject) sym);
3069: }
3070: section->sym = sym;
3071: return(0);
3072: }
3074: /*@
3075: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3077: Not collective
3079: Input Parameters:
3080: . section - the section describing data layout
3082: Output Parameters:
3083: . sym - the symmetry describing the affect of orientation on the access of the data
3085: Level: developer
3087: .seealso: PetscSectionSetSym(), PetscSectionSymCreate()
3088: @*/
3089: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3090: {
3093: *sym = section->sym;
3094: return(0);
3095: }
3097: /*@
3098: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3100: Collective on PetscSection
3102: Input Parameters:
3103: + section - the section describing data layout
3104: . field - the field number
3105: - sym - the symmetry describing the affect of orientation on the access of the data
3107: Level: developer
3109: .seealso: PetscSectionGetFieldSym(), PetscSectionSymCreate()
3110: @*/
3111: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3112: {
3117: 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);
3118: PetscSectionSetSym(section->field[field],sym);
3119: return(0);
3120: }
3122: /*@
3123: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3125: Collective on PetscSection
3127: Input Parameters:
3128: + section - the section describing data layout
3129: - field - the field number
3131: Output Parameters:
3132: . sym - the symmetry describing the affect of orientation on the access of the data
3134: Level: developer
3136: .seealso: PetscSectionSetFieldSym(), PetscSectionSymCreate()
3137: @*/
3138: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3139: {
3142: 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);
3143: *sym = section->field[field]->sym;
3144: return(0);
3145: }
3147: /*@C
3148: PetscSectionGetPointSyms - Get the symmetries for a set of points in a PetscSection under specific orientations.
3150: Not collective
3152: Input Parameters:
3153: + section - the section
3154: . numPoints - the number of points
3155: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3156: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3157: context, see DMPlexGetConeOrientation()).
3159: Output Parameter:
3160: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3161: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3162: identity).
3164: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3165: .vb
3166: const PetscInt **perms;
3167: const PetscScalar **rots;
3168: PetscInt lOffset;
3170: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3171: for (i = 0, lOffset = 0; i < numPoints; i++) {
3172: PetscInt point = points[2*i], dof, sOffset;
3173: const PetscInt *perm = perms ? perms[i] : NULL;
3174: const PetscScalar *rot = rots ? rots[i] : NULL;
3176: PetscSectionGetDof(section,point,&dof);
3177: PetscSectionGetOffset(section,point,&sOffset);
3179: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3180: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3181: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3182: lOffset += dof;
3183: }
3184: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3185: .ve
3187: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3188: .vb
3189: const PetscInt **perms;
3190: const PetscScalar **rots;
3191: PetscInt lOffset;
3193: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3194: for (i = 0, lOffset = 0; i < numPoints; i++) {
3195: PetscInt point = points[2*i], dof, sOffset;
3196: const PetscInt *perm = perms ? perms[i] : NULL;
3197: const PetscScalar *rot = rots ? rots[i] : NULL;
3199: PetscSectionGetDof(section,point,&dof);
3200: PetscSectionGetOffset(section,point,&sOff);
3202: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3203: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3204: offset += dof;
3205: }
3206: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3207: .ve
3209: Level: developer
3211: .seealso: PetscSectionRestorePointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3212: @*/
3213: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3214: {
3215: PetscSectionSym sym;
3216: PetscErrorCode ierr;
3221: if (perms) *perms = NULL;
3222: if (rots) *rots = NULL;
3223: sym = section->sym;
3224: if (sym && (perms || rots)) {
3225: SymWorkLink link;
3227: if (sym->workin) {
3228: link = sym->workin;
3229: sym->workin = sym->workin->next;
3230: } else {
3231: PetscNewLog(sym,&link);
3232: }
3233: if (numPoints > link->numPoints) {
3234: PetscFree2(*(PetscInt***)&link->perms,*(PetscInt***)&link->rots);
3235: PetscMalloc2(numPoints,(PetscInt***)&link->perms,numPoints,(PetscScalar***)&link->rots);
3236: link->numPoints = numPoints;
3237: }
3238: link->next = sym->workout;
3239: sym->workout = link;
3240: PetscArrayzero((PetscInt**)link->perms,numPoints);
3241: PetscArrayzero((PetscInt**)link->rots,numPoints);
3242: (*sym->ops->getpoints) (sym, section, numPoints, points, link->perms, link->rots);
3243: if (perms) *perms = link->perms;
3244: if (rots) *rots = link->rots;
3245: }
3246: return(0);
3247: }
3249: /*@C
3250: PetscSectionRestorePointSyms - Restore the symmetries returned by PetscSectionGetPointSyms()
3252: Not collective
3254: Input Parameters:
3255: + section - the section
3256: . numPoints - the number of points
3257: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3258: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3259: context, see DMPlexGetConeOrientation()).
3261: Output Parameter:
3262: + perms - The permutations for the given orientations: set to NULL at conclusion
3263: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3265: Level: developer
3267: .seealso: PetscSectionGetPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3268: @*/
3269: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3270: {
3271: PetscSectionSym sym;
3275: sym = section->sym;
3276: if (sym && (perms || rots)) {
3277: SymWorkLink *p,link;
3279: for (p=&sym->workout; (link=*p); p=&link->next) {
3280: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3281: *p = link->next;
3282: link->next = sym->workin;
3283: sym->workin = link;
3284: if (perms) *perms = NULL;
3285: if (rots) *rots = NULL;
3286: return(0);
3287: }
3288: }
3289: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
3290: }
3291: return(0);
3292: }
3294: /*@C
3295: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a PetscSection under specific orientations.
3297: Not collective
3299: Input Parameters:
3300: + section - the section
3301: . field - the field of the section
3302: . numPoints - the number of points
3303: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3304: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3305: context, see DMPlexGetConeOrientation()).
3307: Output Parameter:
3308: + perms - The permutations for the given orientations (or NULL if there is no symmetry or the permutation is the identity).
3309: - rots - The field rotations symmetries for the given orientations (or NULL if there is no symmetry or the rotations are all
3310: identity).
3312: Level: developer
3314: .seealso: PetscSectionGetPointSyms(), PetscSectionRestoreFieldPointSyms()
3315: @*/
3316: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3317: {
3322: 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);
3323: PetscSectionGetPointSyms(section->field[field],numPoints,points,perms,rots);
3324: return(0);
3325: }
3327: /*@C
3328: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by PetscSectionGetFieldPointSyms()
3330: Not collective
3332: Input Parameters:
3333: + section - the section
3334: . field - the field number
3335: . numPoints - the number of points
3336: - points - an array of size 2 * numPoints, containing a list of (point, orientation) pairs. (An orientation is an
3337: arbitrary integer: its interpretation is up to sym. Orientations are used by DM: for their interpretation in that
3338: context, see DMPlexGetConeOrientation()).
3340: Output Parameter:
3341: + perms - The permutations for the given orientations: set to NULL at conclusion
3342: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3344: Level: developer
3346: .seealso: PetscSectionRestorePointSyms(), petscSectionGetFieldPointSyms(), PetscSectionSymCreate(), PetscSectionSetSym(), PetscSectionGetSym()
3347: @*/
3348: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3349: {
3354: 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);
3355: PetscSectionRestorePointSyms(section->field[field],numPoints,points,perms,rots);
3356: return(0);
3357: }
3359: /*@
3360: PetscSectionGetUseFieldOffsets - Get the flag to use field offsets directly in a global section, rather than just the point offset
3362: Not collective
3364: Input Parameter:
3365: . s - the global PetscSection
3367: Output Parameters:
3368: . flg - the flag
3370: Level: developer
3372: .seealso: PetscSectionSetChart(), PetscSectionCreate()
3373: @*/
3374: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3375: {
3378: *flg = s->useFieldOff;
3379: return(0);
3380: }
3382: /*@
3383: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3385: Not collective
3387: Input Parameters:
3388: + s - the global PetscSection
3389: - flg - the flag
3391: Level: developer
3393: .seealso: PetscSectionGetUseFieldOffsets(), PetscSectionSetChart(), PetscSectionCreate()
3394: @*/
3395: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3396: {
3399: s->useFieldOff = flg;
3400: return(0);
3401: }
3403: #define PetscSectionExpandPoints_Loop(TYPE) \
3404: { \
3405: PetscInt i, n, o0, o1, size; \
3406: TYPE *a0 = (TYPE*)origArray, *a1; \
3407: PetscSectionGetStorageSize(s, &size); \
3408: PetscMalloc1(size, &a1); \
3409: for (i=0; i<npoints; i++) { \
3410: PetscSectionGetOffset(origSection, points_[i], &o0); \
3411: PetscSectionGetOffset(s, i, &o1); \
3412: PetscSectionGetDof(s, i, &n); \
3413: PetscMemcpy(&a1[o1], &a0[o0], n*unitsize); \
3414: } \
3415: *newArray = (void*)a1; \
3416: }
3418: /*@
3419: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3421: Not collective
3423: Input Parameters:
3424: + origSection - the PetscSection describing the layout of the array
3425: . dataType - MPI_Datatype describing the data type of the array (currently only MPIU_INT, MPIU_SCALAR, MPIU_REAL)
3426: . origArray - the array; its size must be equal to the storage size of origSection
3427: - points - IS with points to extract; its indices must lie in the chart of origSection
3429: Output Parameters:
3430: + newSection - the new PetscSection desribing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3431: - newArray - the array of the extracted DOFs; its size is the storage size of newSection
3433: Level: developer
3435: .seealso: PetscSectionGetChart(), PetscSectionGetDof(), PetscSectionGetStorageSize(), PetscSectionCreate()
3436: @*/
3437: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3438: {
3439: PetscSection s;
3440: const PetscInt *points_;
3441: PetscInt i, n, npoints, pStart, pEnd;
3442: PetscMPIInt unitsize;
3443: PetscErrorCode ierr;
3451: MPI_Type_size(dataType, &unitsize);
3452: ISGetLocalSize(points, &npoints);
3453: ISGetIndices(points, &points_);
3454: PetscSectionGetChart(origSection, &pStart, &pEnd);
3455: PetscSectionCreate(PETSC_COMM_SELF, &s);
3456: PetscSectionSetChart(s, 0, npoints);
3457: for (i=0; i<npoints; i++) {
3458: 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);
3459: PetscSectionGetDof(origSection, points_[i], &n);
3460: PetscSectionSetDof(s, i, n);
3461: }
3462: PetscSectionSetUp(s);
3463: if (newArray) {
3464: if (dataType == MPIU_INT) {PetscSectionExpandPoints_Loop(PetscInt);}
3465: else if (dataType == MPIU_SCALAR) {PetscSectionExpandPoints_Loop(PetscScalar);}
3466: else if (dataType == MPIU_REAL) {PetscSectionExpandPoints_Loop(PetscReal);}
3467: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3468: }
3469: if (newSection) {
3470: *newSection = s;
3471: } else {
3472: PetscSectionDestroy(&s);
3473: }
3474: return(0);
3475: }