Actual source code: section.c
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsf.h>
8: PetscClassId PETSC_SECTION_CLASSID;
10: /*@
11: PetscSectionCreate - Allocates a `PetscSection` and sets the map contents to the default.
13: Collective
15: Input Parameters:
16: + comm - the MPI communicator
17: - s - pointer to the section
19: Level: beginner
21: Notes:
22: Typical calling sequence
23: .vb
24: PetscSectionCreate(MPI_Comm,PetscSection *);!
25: PetscSectionSetNumFields(PetscSection, numFields);
26: PetscSectionSetChart(PetscSection,low,high);
27: PetscSectionSetDof(PetscSection,point,numdof);
28: PetscSectionSetUp(PetscSection);
29: PetscSectionGetOffset(PetscSection,point,PetscInt *);
30: PetscSectionDestroy(PetscSection);
31: .ve
33: The `PetscSection` object and methods are intended to be used in the PETSc `Vec` and `Mat` implementations. The indices returned by the `PetscSection` are appropriate for the kind of `Vec` it is associated with. For example, if the vector being indexed is a local vector, we call the section a local section. If the section indexes a global vector, we call it a global section. For parallel vectors, like global vectors, we use negative indices to indicate dofs owned by other processes.
35: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionDestroy()`, `PetscSectionCreateGlobalSection()`
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
39: PetscFunctionBegin;
40: PetscAssertPointer(s, 2);
41: PetscCall(ISInitializePackage());
43: PetscCall(PetscHeaderCreate(*s, PETSC_SECTION_CLASSID, "PetscSection", "Section", "IS", comm, PetscSectionDestroy, PetscSectionView));
44: (*s)->pStart = -1;
45: (*s)->pEnd = -1;
46: (*s)->perm = NULL;
47: (*s)->pointMajor = PETSC_TRUE;
48: (*s)->includesConstraints = PETSC_TRUE;
49: (*s)->atlasDof = NULL;
50: (*s)->atlasOff = NULL;
51: (*s)->bc = NULL;
52: (*s)->bcIndices = NULL;
53: (*s)->setup = PETSC_FALSE;
54: (*s)->numFields = 0;
55: (*s)->fieldNames = NULL;
56: (*s)->field = NULL;
57: (*s)->useFieldOff = PETSC_FALSE;
58: (*s)->compNames = NULL;
59: (*s)->clObj = NULL;
60: (*s)->clHash = NULL;
61: (*s)->clSection = NULL;
62: (*s)->clPoints = NULL;
63: PetscCall(PetscSectionInvalidateMaxDof_Internal(*s));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*@
68: PetscSectionCopy - Creates a shallow (if possible) copy of the `PetscSection`
70: Collective
72: Input Parameter:
73: . section - the `PetscSection`
75: Output Parameter:
76: . newSection - the copy
78: Level: intermediate
80: Developer Notes:
81: What exactly does shallow mean in this context?
83: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
84: @*/
85: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
86: {
87: PetscFunctionBegin;
90: PetscCall(PetscSectionCopy_Internal(section, newSection, NULL));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: PetscErrorCode PetscSectionCopy_Internal(PetscSection section, PetscSection newSection, PetscBT constrained_dofs)
95: {
96: PetscSectionSym sym;
97: IS perm;
98: PetscInt numFields, f, c, pStart, pEnd, p;
100: PetscFunctionBegin;
103: PetscCall(PetscSectionReset(newSection));
104: PetscCall(PetscSectionGetNumFields(section, &numFields));
105: if (numFields) PetscCall(PetscSectionSetNumFields(newSection, numFields));
106: for (f = 0; f < numFields; ++f) {
107: const char *fieldName = NULL, *compName = NULL;
108: PetscInt numComp = 0;
110: PetscCall(PetscSectionGetFieldName(section, f, &fieldName));
111: PetscCall(PetscSectionSetFieldName(newSection, f, fieldName));
112: PetscCall(PetscSectionGetFieldComponents(section, f, &numComp));
113: PetscCall(PetscSectionSetFieldComponents(newSection, f, numComp));
114: for (c = 0; c < numComp; ++c) {
115: PetscCall(PetscSectionGetComponentName(section, f, c, &compName));
116: PetscCall(PetscSectionSetComponentName(newSection, f, c, compName));
117: }
118: PetscCall(PetscSectionGetFieldSym(section, f, &sym));
119: PetscCall(PetscSectionSetFieldSym(newSection, f, sym));
120: }
121: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
122: PetscCall(PetscSectionSetChart(newSection, pStart, pEnd));
123: PetscCall(PetscSectionGetPermutation(section, &perm));
124: PetscCall(PetscSectionSetPermutation(newSection, perm));
125: PetscCall(PetscSectionGetSym(section, &sym));
126: PetscCall(PetscSectionSetSym(newSection, sym));
127: for (p = pStart; p < pEnd; ++p) {
128: PetscInt dof, cdof, fcdof = 0;
129: PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));
131: PetscCall(PetscSectionGetDof(section, p, &dof));
132: PetscCall(PetscSectionSetDof(newSection, p, dof));
133: if (force_constrained) cdof = dof;
134: else PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
135: if (cdof) PetscCall(PetscSectionSetConstraintDof(newSection, p, cdof));
136: for (f = 0; f < numFields; ++f) {
137: PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
138: PetscCall(PetscSectionSetFieldDof(newSection, p, f, dof));
139: if (cdof) {
140: if (force_constrained) fcdof = dof;
141: else PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
142: if (fcdof) PetscCall(PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof));
143: }
144: }
145: }
146: PetscCall(PetscSectionSetUp(newSection));
147: for (p = pStart; p < pEnd; ++p) {
148: PetscInt off, cdof, fcdof = 0;
149: const PetscInt *cInd;
150: PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));
152: /* Must set offsets in case they do not agree with the prefix sums */
153: PetscCall(PetscSectionGetOffset(section, p, &off));
154: PetscCall(PetscSectionSetOffset(newSection, p, off));
155: PetscCall(PetscSectionGetConstraintDof(newSection, p, &cdof));
156: if (cdof) {
157: if (force_constrained) cInd = NULL;
158: else PetscCall(PetscSectionGetConstraintIndices(section, p, &cInd));
159: PetscCall(PetscSectionSetConstraintIndices(newSection, p, cInd));
160: for (f = 0; f < numFields; ++f) {
161: PetscCall(PetscSectionGetFieldOffset(section, p, f, &off));
162: PetscCall(PetscSectionSetFieldOffset(newSection, p, f, off));
163: PetscCall(PetscSectionGetFieldConstraintDof(newSection, p, f, &fcdof));
164: if (fcdof) {
165: if (force_constrained) cInd = NULL;
166: else PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &cInd));
167: PetscCall(PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd));
168: }
169: }
170: }
171: }
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: /*@
176: PetscSectionClone - Creates a shallow (if possible) copy of the `PetscSection`
178: Collective
180: Input Parameter:
181: . section - the `PetscSection`
183: Output Parameter:
184: . newSection - the copy
186: Level: beginner
188: Developer Notes:
189: With standard PETSc terminology this should be called `PetscSectionDuplicate()`
191: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
192: @*/
193: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
194: {
195: PetscFunctionBegin;
197: PetscAssertPointer(newSection, 2);
198: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection));
199: PetscCall(PetscSectionCopy(section, *newSection));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: PetscSectionSetFromOptions - sets parameters in a `PetscSection` from the options database
206: Collective
208: Input Parameter:
209: . s - the `PetscSection`
211: Options Database Key:
212: . -petscsection_point_major - `PETSC_TRUE` for point-major order
214: Level: intermediate
216: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
217: @*/
218: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
219: {
220: PetscFunctionBegin;
222: PetscObjectOptionsBegin((PetscObject)s);
223: PetscCall(PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL));
224: /* process any options handlers added with PetscObjectAddOptionsHandler() */
225: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject));
226: PetscOptionsEnd();
227: PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view"));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: PetscSectionCompare - Compares two sections
234: Collective
236: Input Parameters:
237: + s1 - the first `PetscSection`
238: - s2 - the second `PetscSection`
240: Output Parameter:
241: . congruent - `PETSC_TRUE` if the two sections are congruent, `PETSC_FALSE` otherwise
243: Level: intermediate
245: Note:
246: Field names are disregarded.
248: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
249: @*/
250: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
251: {
252: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
253: const PetscInt *idx1, *idx2;
254: IS perm1, perm2;
255: PetscMPIInt mflg;
257: PetscFunctionBegin;
260: PetscAssertPointer(congruent, 3);
261: *congruent = PETSC_FALSE;
263: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg));
264: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
265: *congruent = PETSC_FALSE;
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: PetscCall(PetscSectionGetChart(s1, &pStart, &pEnd));
270: PetscCall(PetscSectionGetChart(s2, &n1, &n2));
271: if (pStart != n1 || pEnd != n2) goto not_congruent;
273: PetscCall(PetscSectionGetPermutation(s1, &perm1));
274: PetscCall(PetscSectionGetPermutation(s2, &perm2));
275: if (perm1 && perm2) {
276: PetscCall(ISEqual(perm1, perm2, congruent));
277: if (!*congruent) goto not_congruent;
278: } else if (perm1 != perm2) goto not_congruent;
280: for (p = pStart; p < pEnd; ++p) {
281: PetscCall(PetscSectionGetOffset(s1, p, &n1));
282: PetscCall(PetscSectionGetOffset(s2, p, &n2));
283: if (n1 != n2) goto not_congruent;
285: PetscCall(PetscSectionGetDof(s1, p, &n1));
286: PetscCall(PetscSectionGetDof(s2, p, &n2));
287: if (n1 != n2) goto not_congruent;
289: PetscCall(PetscSectionGetConstraintDof(s1, p, &ncdof));
290: PetscCall(PetscSectionGetConstraintDof(s2, p, &n2));
291: if (ncdof != n2) goto not_congruent;
293: PetscCall(PetscSectionGetConstraintIndices(s1, p, &idx1));
294: PetscCall(PetscSectionGetConstraintIndices(s2, p, &idx2));
295: PetscCall(PetscArraycmp(idx1, idx2, ncdof, congruent));
296: if (!*congruent) goto not_congruent;
297: }
299: PetscCall(PetscSectionGetNumFields(s1, &nfields));
300: PetscCall(PetscSectionGetNumFields(s2, &n2));
301: if (nfields != n2) goto not_congruent;
303: for (f = 0; f < nfields; ++f) {
304: PetscCall(PetscSectionGetFieldComponents(s1, f, &n1));
305: PetscCall(PetscSectionGetFieldComponents(s2, f, &n2));
306: if (n1 != n2) goto not_congruent;
308: for (p = pStart; p < pEnd; ++p) {
309: PetscCall(PetscSectionGetFieldOffset(s1, p, f, &n1));
310: PetscCall(PetscSectionGetFieldOffset(s2, p, f, &n2));
311: if (n1 != n2) goto not_congruent;
313: PetscCall(PetscSectionGetFieldDof(s1, p, f, &n1));
314: PetscCall(PetscSectionGetFieldDof(s2, p, f, &n2));
315: if (n1 != n2) goto not_congruent;
317: PetscCall(PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof));
318: PetscCall(PetscSectionGetFieldConstraintDof(s2, p, f, &n2));
319: if (nfcdof != n2) goto not_congruent;
321: PetscCall(PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1));
322: PetscCall(PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2));
323: PetscCall(PetscArraycmp(idx1, idx2, nfcdof, congruent));
324: if (!*congruent) goto not_congruent;
325: }
326: }
328: *congruent = PETSC_TRUE;
329: not_congruent:
330: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, congruent, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1)));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: /*@
335: PetscSectionGetNumFields - Returns the number of fields in a `PetscSection`, or 0 if no fields were defined.
337: Not Collective
339: Input Parameter:
340: . s - the `PetscSection`
342: Output Parameter:
343: . numFields - the number of fields defined, or 0 if none were defined
345: Level: intermediate
347: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
348: @*/
349: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
350: {
351: PetscFunctionBegin;
353: PetscAssertPointer(numFields, 2);
354: *numFields = s->numFields;
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: /*@
359: PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`
361: Not Collective
363: Input Parameters:
364: + s - the `PetscSection`
365: - numFields - the number of fields
367: Level: intermediate
369: Notes:
370: Calling this destroys all the information in the `PetscSection` including the chart.
372: You must call `PetscSectionSetChart()` after calling this.
374: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetNumFields()`, `PetscSectionSetChart()`, `PetscSectionReset()`
375: @*/
376: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
377: {
378: PetscFunctionBegin;
380: PetscCheck(numFields > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %" PetscInt_FMT " must be positive", numFields);
381: PetscCall(PetscSectionReset(s));
383: s->numFields = numFields;
384: PetscCall(PetscMalloc1(s->numFields, &s->numFieldComponents));
385: PetscCall(PetscMalloc1(s->numFields, &s->fieldNames));
386: PetscCall(PetscMalloc1(s->numFields, &s->compNames));
387: PetscCall(PetscMalloc1(s->numFields, &s->field));
388: for (PetscInt f = 0; f < s->numFields; ++f) {
389: char name[64];
391: s->numFieldComponents[f] = 1;
393: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]));
394: PetscCall(PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f));
395: PetscCall(PetscStrallocpy(name, &s->fieldNames[f]));
396: PetscCall(PetscSNPrintf(name, 64, "Component_0"));
397: PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]));
398: PetscCall(PetscStrallocpy(name, &s->compNames[f][0]));
399: }
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: /*@
404: PetscSectionGetFieldName - Returns the name of a field in the `PetscSection`
406: Not Collective
408: Input Parameters:
409: + s - the `PetscSection`
410: - field - the field number
412: Output Parameter:
413: . fieldName - the field name
415: Level: intermediate
417: Note:
418: Will error if the field number is out of range
420: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
421: @*/
422: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
423: {
424: PetscFunctionBegin;
426: PetscAssertPointer(fieldName, 3);
427: PetscSectionCheckValidField(field, s->numFields);
428: *fieldName = s->fieldNames[field];
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
432: /*@
433: PetscSectionSetFieldName - Sets the name of a field in the `PetscSection`
435: Not Collective
437: Input Parameters:
438: + s - the `PetscSection`
439: . field - the field number
440: - fieldName - the field name
442: Level: intermediate
444: Note:
445: Will error if the field number is out of range
447: .seealso: [PetscSection](ch_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
448: @*/
449: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
450: {
451: PetscFunctionBegin;
453: if (fieldName) PetscAssertPointer(fieldName, 3);
454: PetscSectionCheckValidField(field, s->numFields);
455: PetscCall(PetscFree(s->fieldNames[field]));
456: PetscCall(PetscStrallocpy(fieldName, &s->fieldNames[field]));
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /*@
461: PetscSectionGetComponentName - Gets the name of a field component in the `PetscSection`
463: Not Collective
465: Input Parameters:
466: + s - the `PetscSection`
467: . field - the field number
468: - comp - the component number
470: Output Parameter:
471: . compName - the component name
473: Level: intermediate
475: Note:
476: Will error if the field or component number do not exist
478: Developer Notes:
479: The function name should have Field in it since they are field components.
481: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
482: `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
483: @*/
484: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
485: {
486: PetscFunctionBegin;
488: PetscAssertPointer(compName, 4);
489: PetscSectionCheckValidField(field, s->numFields);
490: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
491: *compName = s->compNames[field][comp];
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@
496: PetscSectionSetComponentName - Sets the name of a field component in the `PetscSection`
498: Not Collective
500: Input Parameters:
501: + s - the `PetscSection`
502: . field - the field number
503: . comp - the component number
504: - compName - the component name
506: Level: advanced
508: Note:
509: Will error if the field or component number do not exist
511: Developer Notes:
512: The function name should have Field in it since they are field components.
514: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
515: `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
516: @*/
517: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
518: {
519: PetscFunctionBegin;
521: if (compName) PetscAssertPointer(compName, 4);
522: PetscSectionCheckValidField(field, s->numFields);
523: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
524: PetscCall(PetscFree(s->compNames[field][comp]));
525: PetscCall(PetscStrallocpy(compName, &s->compNames[field][comp]));
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
532: Not Collective
534: Input Parameters:
535: + s - the `PetscSection`
536: - field - the field number
538: Output Parameter:
539: . numComp - the number of field components
541: Level: advanced
543: Developer Notes:
544: This function is misnamed. There is a Num in `PetscSectionGetNumFields()` but not in this name
546: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`,
547: `PetscSectionSetComponentName()`, `PetscSectionGetComponentName()`
548: @*/
549: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
550: {
551: PetscFunctionBegin;
553: PetscAssertPointer(numComp, 3);
554: PetscSectionCheckValidField(field, s->numFields);
555: *numComp = s->numFieldComponents[field];
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: /*@
560: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
562: Not Collective
564: Input Parameters:
565: + s - the `PetscSection`
566: . field - the field number
567: - numComp - the number of field components
569: Level: advanced
571: Note:
572: This number can be different than the values set with `PetscSectionSetFieldDof()`. It can be used to indicate the number of
573: components in the field of the underlying physical model which may be different than the number of degrees of freedom needed
574: at a point in a discretization. For example, if in three dimensions the field is velocity, it will have 3 components, u, v, and w but
575: an face based model for velocity (where the velocity normal to the face is stored) there is only 1 dof for each face point.
577: The value set with this function are not needed or used in `PetscSectionSetUp()`.
579: Developer Notes:
580: This function is misnamed. There is a Num in `PetscSectionSetNumFields()` but not in this name
582: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionSetComponentName()`,
583: `PetscSectionGetComponentName()`, `PetscSectionGetNumFields()`
584: @*/
585: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
586: {
587: PetscFunctionBegin;
589: PetscSectionCheckValidField(field, s->numFields);
590: if (s->compNames) {
591: for (PetscInt c = 0; c < s->numFieldComponents[field]; ++c) PetscCall(PetscFree(s->compNames[field][c]));
592: PetscCall(PetscFree(s->compNames[field]));
593: }
595: s->numFieldComponents[field] = numComp;
596: if (numComp) {
597: PetscCall(PetscMalloc1(numComp, &s->compNames[field]));
598: for (PetscInt c = 0; c < numComp; ++c) {
599: char name[64];
601: PetscCall(PetscSNPrintf(name, 64, "%" PetscInt_FMT, c));
602: PetscCall(PetscStrallocpy(name, &s->compNames[field][c]));
603: }
604: }
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: /*@
609: PetscSectionGetChart - Returns the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
611: Not Collective
613: Input Parameter:
614: . s - the `PetscSection`
616: Output Parameters:
617: + pStart - the first point
618: - pEnd - one past the last point
620: Level: intermediate
622: Note:
623: The chart may be thought of as the bounds on the points (indices) one may use to index into numerical data that is associated with
624: the `PetscSection` data layout.
626: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
627: @*/
628: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
629: {
630: PetscFunctionBegin;
632: if (pStart) *pStart = s->pStart;
633: if (pEnd) *pEnd = s->pEnd;
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: /*@
638: PetscSectionSetChart - Sets the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
640: Not Collective
642: Input Parameters:
643: + s - the `PetscSection`
644: . pStart - the first `point`
645: - pEnd - one past the last point, `pStart` $ \le $ `pEnd`
647: Level: intermediate
649: Notes:
650: The chart may be thought of as the bounds on the points (indices) one may use to index into numerical data that is associated with
651: the `PetscSection` data layout.
653: The charts on different MPI processes may (and often do) overlap
655: If you intend to use `PetscSectionSetNumFields()` it must be called before this call.
657: The chart for all fields created with `PetscSectionSetNumFields()` is the same as this chart.
659: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`, `PetscSectionSetNumFields()`
660: @*/
661: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
662: {
663: PetscFunctionBegin;
665: PetscCheck(pEnd >= pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Chart pEnd %" PetscInt_FMT " cannot be smaller than chart pStart %" PetscInt_FMT, pEnd, pStart);
666: if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(PETSC_SUCCESS);
667: /* Cannot Reset() because it destroys field information */
668: s->setup = PETSC_FALSE;
669: PetscCall(PetscSectionDestroy(&s->bc));
670: PetscCall(PetscFree(s->bcIndices));
671: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
673: s->pStart = pStart;
674: s->pEnd = pEnd;
675: PetscCall(PetscMalloc2(pEnd - pStart, &s->atlasDof, pEnd - pStart, &s->atlasOff));
676: PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart));
677: for (PetscInt f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: /*@
682: PetscSectionGetPermutation - Returns the permutation of [0, `pEnd` - `pStart`) or `NULL` that was set with `PetscSectionSetPermutation()`
684: Not Collective
686: Input Parameter:
687: . s - the `PetscSection`
689: Output Parameter:
690: . perm - The permutation as an `IS`
692: Level: intermediate
694: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
695: @*/
696: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
697: {
698: PetscFunctionBegin;
700: if (perm) {
701: PetscAssertPointer(perm, 2);
702: *perm = s->perm;
703: }
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: /*@
708: PetscSectionSetPermutation - Sets a permutation of the chart for this section, [0, `pEnd` - `pStart`), which determines the order to store the `PetscSection` information
710: Not Collective
712: Input Parameters:
713: + s - the `PetscSection`
714: - perm - the permutation of points
716: Level: intermediate
718: Notes:
719: The permutation must be provided before `PetscSectionSetUp()`.
721: The data in the `PetscSection` are permuted but the access via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` is not changed
723: Compare to `PetscSectionPermute()`
725: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetUp()`, `PetscSectionGetPermutation()`, `PetscSectionPermute()`, `PetscSectionCreate()`
726: @*/
727: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
728: {
729: PetscFunctionBegin;
732: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
733: if (s->perm != perm) {
734: PetscCall(ISDestroy(&s->perm));
735: if (perm) {
736: s->perm = perm;
737: PetscCall(PetscObjectReference((PetscObject)s->perm));
738: }
739: }
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: PetscSectionGetBlockStarts - Returns a table indicating which points start new blocks
746: Not Collective, No Fortran Support
748: Input Parameter:
749: . s - the `PetscSection`
751: Output Parameter:
752: . blockStarts - The `PetscBT` with a 1 for each point that begins a block
754: Notes:
755: The table is on [0, `pEnd` - `pStart`).
757: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
759: Level: intermediate
761: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
762: @*/
763: PetscErrorCode PetscSectionGetBlockStarts(PetscSection s, PetscBT *blockStarts)
764: {
765: PetscFunctionBegin;
767: if (blockStarts) {
768: PetscAssertPointer(blockStarts, 2);
769: *blockStarts = s->blockStarts;
770: }
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: /*@
775: PetscSectionSetBlockStarts - Sets a table indicating which points start new blocks
777: Not Collective, No Fortran Support
779: Input Parameters:
780: + s - the `PetscSection`
781: - blockStarts - The `PetscBT` with a 1 for each point that begins a block
783: Level: intermediate
785: Notes:
786: The table is on [0, `pEnd` - `pStart`). PETSc takes ownership of the `PetscBT` when it is passed in and will destroy it. The user should not destroy it.
788: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
790: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionGetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
791: @*/
792: PetscErrorCode PetscSectionSetBlockStarts(PetscSection s, PetscBT blockStarts)
793: {
794: PetscFunctionBegin;
796: if (s->blockStarts != blockStarts) {
797: PetscCall(PetscBTDestroy(&s->blockStarts));
798: s->blockStarts = blockStarts;
799: }
800: PetscFunctionReturn(PETSC_SUCCESS);
801: }
803: /*@
804: PetscSectionGetPointMajor - Returns the flag for dof ordering, `PETSC_TRUE` if it is point major, `PETSC_FALSE` if it is field major
806: Not Collective
808: Input Parameter:
809: . s - the `PetscSection`
811: Output Parameter:
812: . pm - the flag for point major ordering
814: Level: intermediate
816: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetPointMajor()`
817: @*/
818: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
819: {
820: PetscFunctionBegin;
822: PetscAssertPointer(pm, 2);
823: *pm = s->pointMajor;
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: /*@
828: PetscSectionSetPointMajor - Sets the flag for dof ordering, `PETSC_TRUE` for point major, otherwise it will be field major
830: Not Collective
832: Input Parameters:
833: + s - the `PetscSection`
834: - pm - the flag for point major ordering
836: Level: intermediate
838: Note:
839: Field-major order is not recommended unless you are managing the entire problem yourself, since many higher-level functions in PETSc depend on point-major order.
841: Point major order means the degrees of freedom are stored as follows
842: .vb
843: all the degrees of freedom for each point are stored contiguously, one point after another (respecting a permutation set with `PetscSectionSetPermutation()`)
844: for each point
845: the degrees of freedom for each field (starting with the unnamed default field) are listed in order by field
846: .ve
848: Field major order means the degrees of freedom are stored as follows
849: .vb
850: all degrees of freedom for each field (including the unnamed default field) are stored contiguously, one field after another
851: for each field (started with unnamed default field)
852: the degrees of freedom for each point are listed in order by point (respecting a permutation set with `PetscSectionSetPermutation()`)
853: .ve
855: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`, `PetscSectionSetPermutation()`
856: @*/
857: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
858: {
859: PetscFunctionBegin;
861: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
862: s->pointMajor = pm;
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: /*@
867: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets in the `PetscSection`.
868: The value is set with `PetscSectionSetIncludesConstraints()`
870: Not Collective
872: Input Parameter:
873: . s - the `PetscSection`
875: Output Parameter:
876: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
878: Level: intermediate
880: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
881: @*/
882: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
883: {
884: PetscFunctionBegin;
886: PetscAssertPointer(includesConstraints, 2);
887: *includesConstraints = s->includesConstraints;
888: PetscFunctionReturn(PETSC_SUCCESS);
889: }
891: /*@
892: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
894: Not Collective
896: Input Parameters:
897: + s - the `PetscSection`
898: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
900: Level: intermediate
902: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
903: @*/
904: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
905: {
906: PetscFunctionBegin;
908: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
909: s->includesConstraints = includesConstraints;
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: /*@
914: PetscSectionGetDof - Return the total number of degrees of freedom associated with a given point.
916: Not Collective
918: Input Parameters:
919: + s - the `PetscSection`
920: - point - the point
922: Output Parameter:
923: . numDof - the number of dof
925: Level: intermediate
927: Notes:
928: In a global section, this size will be negative for points not owned by this process.
930: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
932: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
933: @*/
934: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
935: {
936: PetscFunctionBeginHot;
938: PetscAssertPointer(numDof, 3);
939: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
940: *numDof = s->atlasDof[point - s->pStart];
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: /*@
945: PetscSectionSetDof - Sets the total number of degrees of freedom associated with a given point.
947: Not Collective
949: Input Parameters:
950: + s - the `PetscSection`
951: . point - the point
952: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
954: Level: intermediate
956: Note:
957: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
959: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
960: @*/
961: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
962: {
963: PetscFunctionBegin;
965: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
966: s->atlasDof[point - s->pStart] = numDof;
967: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@
972: PetscSectionAddDof - Adds to the total number of degrees of freedom associated with a given point.
974: Not Collective
976: Input Parameters:
977: + s - the `PetscSection`
978: . point - the point
979: - numDof - the number of additional dof
981: Level: intermediate
983: Note:
984: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
986: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
987: @*/
988: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
989: {
990: PetscFunctionBeginHot;
992: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
993: PetscCheck(numDof >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numDof %" PetscInt_FMT " should not be negative", numDof);
994: s->atlasDof[point - s->pStart] += numDof;
995: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: /*@
1000: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
1002: Not Collective
1004: Input Parameters:
1005: + s - the `PetscSection`
1006: . point - the point
1007: - field - the field
1009: Output Parameter:
1010: . numDof - the number of dof
1012: Level: intermediate
1014: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
1015: @*/
1016: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1017: {
1018: PetscFunctionBegin;
1020: PetscAssertPointer(numDof, 4);
1021: PetscSectionCheckValidField(field, s->numFields);
1022: PetscCall(PetscSectionGetDof(s->field[field], point, numDof));
1023: PetscFunctionReturn(PETSC_SUCCESS);
1024: }
1026: /*@
1027: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
1029: Not Collective
1031: Input Parameters:
1032: + s - the `PetscSection`
1033: . point - the point
1034: . field - the field
1035: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
1037: Level: intermediate
1039: Note:
1040: When setting the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
1041: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1043: This is equivalent to
1044: .vb
1045: PetscSection fs;
1046: PetscSectionGetField(s,field,&fs)
1047: PetscSectionSetDof(fs,numDof)
1048: .ve
1050: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`, `PetscSectionAddDof()`, `PetscSectionSetDof()`
1051: @*/
1052: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1053: {
1054: PetscFunctionBegin;
1056: PetscSectionCheckValidField(field, s->numFields);
1057: PetscCall(PetscSectionSetDof(s->field[field], point, numDof));
1058: PetscFunctionReturn(PETSC_SUCCESS);
1059: }
1061: /*@
1062: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
1064: Not Collective
1066: Input Parameters:
1067: + s - the `PetscSection`
1068: . point - the point
1069: . field - the field
1070: - numDof - the number of dof
1072: Level: intermediate
1074: Notes:
1075: When adding to the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
1076: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1078: This is equivalent to
1079: .vb
1080: PetscSection fs;
1081: PetscSectionGetField(s,field,&fs)
1082: PetscSectionAddDof(fs,numDof)
1083: .ve
1085: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
1086: @*/
1087: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1088: {
1089: PetscFunctionBegin;
1091: PetscSectionCheckValidField(field, s->numFields);
1092: PetscCall(PetscSectionAddDof(s->field[field], point, numDof));
1093: PetscFunctionReturn(PETSC_SUCCESS);
1094: }
1096: /*@
1097: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
1099: Not Collective
1101: Input Parameters:
1102: + s - the `PetscSection`
1103: - point - the point
1105: Output Parameter:
1106: . numDof - the number of dof which are fixed by constraints
1108: Level: intermediate
1110: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
1111: @*/
1112: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
1113: {
1114: PetscFunctionBegin;
1116: PetscAssertPointer(numDof, 3);
1117: if (s->bc) PetscCall(PetscSectionGetDof(s->bc, point, numDof));
1118: else *numDof = 0;
1119: PetscFunctionReturn(PETSC_SUCCESS);
1120: }
1122: /*@
1123: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
1125: Not Collective
1127: Input Parameters:
1128: + s - the `PetscSection`
1129: . point - the point
1130: - numDof - the number of dof which are fixed by constraints
1132: Level: intermediate
1134: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1135: @*/
1136: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1137: {
1138: PetscFunctionBegin;
1140: if (numDof) {
1141: PetscCall(PetscSectionCheckConstraints_Private(s));
1142: PetscCall(PetscSectionSetDof(s->bc, point, numDof));
1143: }
1144: PetscFunctionReturn(PETSC_SUCCESS);
1145: }
1147: /*@
1148: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
1150: Not Collective
1152: Input Parameters:
1153: + s - the `PetscSection`
1154: . point - the point
1155: - numDof - the number of additional dof which are fixed by constraints
1157: Level: intermediate
1159: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1160: @*/
1161: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1162: {
1163: PetscFunctionBegin;
1165: if (numDof) {
1166: PetscCall(PetscSectionCheckConstraints_Private(s));
1167: PetscCall(PetscSectionAddDof(s->bc, point, numDof));
1168: }
1169: PetscFunctionReturn(PETSC_SUCCESS);
1170: }
1172: /*@
1173: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
1175: Not Collective
1177: Input Parameters:
1178: + s - the `PetscSection`
1179: . point - the point
1180: - field - the field
1182: Output Parameter:
1183: . numDof - the number of dof which are fixed by constraints
1185: Level: intermediate
1187: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1188: @*/
1189: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1190: {
1191: PetscFunctionBegin;
1193: PetscAssertPointer(numDof, 4);
1194: PetscSectionCheckValidField(field, s->numFields);
1195: PetscCall(PetscSectionGetConstraintDof(s->field[field], point, numDof));
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: /*@
1200: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1202: Not Collective
1204: Input Parameters:
1205: + s - the `PetscSection`
1206: . point - the point
1207: . field - the field
1208: - numDof - the number of dof which are fixed by constraints
1210: Level: intermediate
1212: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1213: @*/
1214: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1215: {
1216: PetscFunctionBegin;
1218: PetscSectionCheckValidField(field, s->numFields);
1219: PetscCall(PetscSectionSetConstraintDof(s->field[field], point, numDof));
1220: PetscFunctionReturn(PETSC_SUCCESS);
1221: }
1223: /*@
1224: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1226: Not Collective
1228: Input Parameters:
1229: + s - the `PetscSection`
1230: . point - the point
1231: . field - the field
1232: - numDof - the number of additional dof which are fixed by constraints
1234: Level: intermediate
1236: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1237: @*/
1238: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1239: {
1240: PetscFunctionBegin;
1242: PetscSectionCheckValidField(field, s->numFields);
1243: PetscCall(PetscSectionAddConstraintDof(s->field[field], point, numDof));
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: /*@
1248: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1250: Not Collective
1252: Input Parameter:
1253: . s - the `PetscSection`
1255: Level: advanced
1257: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1258: @*/
1259: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1260: {
1261: PetscFunctionBegin;
1263: if (s->bc) {
1264: const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;
1266: PetscCall(PetscSectionSetUp(s->bc));
1267: if (last >= 0) PetscCall(PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices));
1268: else s->bcIndices = NULL;
1269: }
1270: PetscFunctionReturn(PETSC_SUCCESS);
1271: }
1273: /*@
1274: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point in preparation for use of the `PetscSection`
1276: Not Collective
1278: Input Parameter:
1279: . s - the `PetscSection`
1281: Level: intermediate
1283: Notes:
1284: If used, `PetscSectionSetPermutation()` must be called before this routine.
1286: `PetscSectionSetPointMajor()`, cannot be called after this routine.
1288: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
1289: @*/
1290: PetscErrorCode PetscSectionSetUp(PetscSection s)
1291: {
1292: PetscInt f;
1293: const PetscInt *pind = NULL;
1294: PetscCount offset = 0;
1296: PetscFunctionBegin;
1298: if (s->setup) PetscFunctionReturn(PETSC_SUCCESS);
1299: s->setup = PETSC_TRUE;
1300: /* Set offsets and field offsets for all points */
1301: /* Assume that all fields have the same chart */
1302: PetscCheck(s->includesConstraints, PETSC_COMM_SELF, PETSC_ERR_SUP, "PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE");
1303: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1304: if (s->pointMajor) {
1305: PetscCount foff;
1306: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1307: const PetscInt q = pind ? pind[p] : p;
1309: /* Set point offset */
1310: PetscCall(PetscIntCast(offset, &s->atlasOff[q]));
1311: offset += s->atlasDof[q];
1312: /* Set field offset */
1313: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1314: PetscSection sf = s->field[f];
1316: PetscCall(PetscIntCast(foff, &sf->atlasOff[q]));
1317: foff += sf->atlasDof[q];
1318: }
1319: }
1320: } else {
1321: /* Set field offsets for all points */
1322: for (f = 0; f < s->numFields; ++f) {
1323: PetscSection sf = s->field[f];
1325: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1326: const PetscInt q = pind ? pind[p] : p;
1328: PetscCall(PetscIntCast(offset, &sf->atlasOff[q]));
1329: offset += sf->atlasDof[q];
1330: }
1331: }
1332: /* Disable point offsets since these are unused */
1333: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1334: }
1335: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1336: /* Setup BC sections */
1337: PetscCall(PetscSectionSetUpBC(s));
1338: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetUpBC(s->field[f]));
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: /*@
1343: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the `PetscSection`
1345: Not Collective
1347: Input Parameter:
1348: . s - the `PetscSection`
1350: Output Parameter:
1351: . maxDof - the maximum dof
1353: Level: intermediate
1355: Notes:
1356: The returned number is up-to-date without need for `PetscSectionSetUp()`.
1358: This is the maximum over all points of the sum of the number of dof in the unnamed default field plus all named fields. This is equivalent to
1359: the maximum over all points of the value returned by `PetscSectionGetDof()` on this MPI process
1361: Developer Notes:
1362: The returned number is calculated lazily and stashed.
1364: A call to `PetscSectionInvalidateMaxDof_Internal()` invalidates the stashed value.
1366: `PetscSectionInvalidateMaxDof_Internal()` is called in `PetscSectionSetDof()`, `PetscSectionAddDof()` and `PetscSectionReset()`
1368: It should also be called every time `atlasDof` is modified directly.
1370: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1371: @*/
1372: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1373: {
1374: PetscInt p;
1376: PetscFunctionBegin;
1378: PetscAssertPointer(maxDof, 2);
1379: if (s->maxDof == PETSC_INT_MIN) {
1380: s->maxDof = 0;
1381: for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1382: }
1383: *maxDof = s->maxDof;
1384: PetscFunctionReturn(PETSC_SUCCESS);
1385: }
1387: /*@
1388: PetscSectionGetStorageSize - Return the size of an array or local `Vec` capable of holding all the degrees of freedom defined in a `PetscSection`
1390: Not Collective
1392: Input Parameter:
1393: . s - the `PetscSection`
1395: Output Parameter:
1396: . size - the size of an array which can hold all the dofs
1398: Level: intermediate
1400: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1401: @*/
1402: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1403: {
1404: PetscInt64 n = 0;
1406: PetscFunctionBegin;
1408: PetscAssertPointer(size, 2);
1409: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1410: PetscCall(PetscIntCast(n, size));
1411: PetscFunctionReturn(PETSC_SUCCESS);
1412: }
1414: /*@
1415: PetscSectionGetConstrainedStorageSize - Return the size of an array or local `Vec` capable of holding all unconstrained degrees of freedom in a `PetscSection`
1417: Not Collective
1419: Input Parameter:
1420: . s - the `PetscSection`
1422: Output Parameter:
1423: . size - the size of an array which can hold all unconstrained dofs
1425: Level: intermediate
1427: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1428: @*/
1429: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1430: {
1431: PetscInt64 n = 0;
1433: PetscFunctionBegin;
1435: PetscAssertPointer(size, 2);
1436: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1437: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1438: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1439: }
1440: PetscCall(PetscIntCast(n, size));
1441: PetscFunctionReturn(PETSC_SUCCESS);
1442: }
1444: /*@
1445: PetscSectionCreateGlobalSection - Create a parallel section describing the global layout using
1446: a local (sequential) `PetscSection` on each MPI process and a `PetscSF` describing the section point overlap.
1448: Input Parameters:
1449: + s - The `PetscSection` for the local field layout
1450: . sf - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1451: . usePermutation - By default this is `PETSC_TRUE`, meaning any permutation of the local section is transferred to the global section
1452: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1453: - localOffsets - If `PETSC_TRUE`, use local rather than global offsets for the points
1455: Output Parameter:
1456: . gsection - The `PetscSection` for the global field layout
1458: Level: intermediate
1460: Notes:
1461: On each MPI process `gsection` inherits the chart of the `s` on that process.
1463: This sets negative sizes and offsets to points not owned by this process as defined by `sf` but that are within the local value of the chart of `gsection`.
1464: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1466: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCreateGlobalSectionCensored()`
1467: @*/
1468: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool usePermutation, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1469: {
1470: PetscSection gs;
1471: const PetscInt *pind = NULL;
1472: PetscInt *recv = NULL, *neg = NULL;
1473: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1474: PetscInt numFields, f, numComponents;
1475: PetscInt foff;
1477: PetscFunctionBegin;
1483: PetscAssertPointer(gsection, 6);
1484: PetscCheck(s->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
1485: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs));
1486: PetscCall(PetscSectionGetNumFields(s, &numFields));
1487: if (numFields > 0) PetscCall(PetscSectionSetNumFields(gs, numFields));
1488: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1489: PetscCall(PetscSectionSetChart(gs, pStart, pEnd));
1490: gs->includesConstraints = includeConstraints;
1491: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1492: nlocal = nroots; /* The local/leaf space matches global/root space */
1493: /* Must allocate for all points visible to SF, which may be more than this section */
1494: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1495: PetscCall(PetscSFGetLeafRange(sf, NULL, &maxleaf));
1496: PetscCheck(nroots >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %" PetscInt_FMT " < pEnd %" PetscInt_FMT, nroots, pEnd);
1497: PetscCheck(maxleaf < nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %" PetscInt_FMT " >= nroots %" PetscInt_FMT, maxleaf, nroots);
1498: PetscCall(PetscMalloc2(nroots, &neg, nlocal, &recv));
1499: PetscCall(PetscArrayzero(neg, nroots));
1500: }
1501: /* Mark all local points with negative dof */
1502: for (p = pStart; p < pEnd; ++p) {
1503: PetscCall(PetscSectionGetDof(s, p, &dof));
1504: PetscCall(PetscSectionSetDof(gs, p, dof));
1505: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1506: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(gs, p, cdof));
1507: if (neg) neg[p] = -(dof + 1);
1508: }
1509: PetscCall(PetscSectionSetUpBC(gs));
1510: if (gs->bcIndices) PetscCall(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]));
1511: if (nroots >= 0) {
1512: PetscCall(PetscArrayzero(recv, nlocal));
1513: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1514: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1515: for (p = pStart; p < pEnd; ++p) {
1516: if (recv[p] < 0) {
1517: gs->atlasDof[p - pStart] = recv[p];
1518: PetscCall(PetscSectionGetDof(s, p, &dof));
1519: PetscCheck(-(recv[p] + 1) == dof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global dof %" PetscInt_FMT " for point %" PetscInt_FMT " is not the unconstrained %" PetscInt_FMT, -(recv[p] + 1), p, dof);
1520: }
1521: }
1522: }
1523: /* Calculate new sizes, get process offset, and calculate point offsets */
1524: if (usePermutation && s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1525: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1526: const PetscInt q = pind ? pind[p] : p;
1528: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1529: gs->atlasOff[q] = off;
1530: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1531: }
1532: if (!localOffsets) {
1533: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)sf)));
1534: globalOff -= off;
1535: }
1536: for (p = pStart, off = 0; p < pEnd; ++p) {
1537: gs->atlasOff[p - pStart] += globalOff;
1538: if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1539: }
1540: if (usePermutation && s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1541: /* Put in negative offsets for ghost points */
1542: if (nroots >= 0) {
1543: PetscCall(PetscArrayzero(recv, nlocal));
1544: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1545: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1546: for (p = pStart; p < pEnd; ++p) {
1547: if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1548: }
1549: }
1550: PetscCall(PetscFree2(neg, recv));
1551: /* Set field dofs/offsets/constraints */
1552: for (f = 0; f < numFields; ++f) {
1553: const char *name;
1555: gs->field[f]->includesConstraints = includeConstraints;
1556: PetscCall(PetscSectionGetFieldComponents(s, f, &numComponents));
1557: PetscCall(PetscSectionSetFieldComponents(gs, f, numComponents));
1558: PetscCall(PetscSectionGetFieldName(s, f, &name));
1559: PetscCall(PetscSectionSetFieldName(gs, f, name));
1560: }
1561: for (p = pStart; p < pEnd; ++p) {
1562: PetscCall(PetscSectionGetOffset(gs, p, &off));
1563: for (f = 0, foff = off; f < numFields; ++f) {
1564: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
1565: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetFieldConstraintDof(gs, p, f, cdof));
1566: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
1567: PetscCall(PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof));
1568: PetscCall(PetscSectionSetFieldOffset(gs, p, f, foff));
1569: PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof));
1570: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1571: }
1572: }
1573: for (f = 0; f < numFields; ++f) {
1574: PetscSection gfs = gs->field[f];
1576: PetscCall(PetscSectionSetUpBC(gfs));
1577: if (gfs->bcIndices) PetscCall(PetscArraycpy(gfs->bcIndices, s->field[f]->bcIndices, gfs->bc->atlasOff[gfs->bc->pEnd - gfs->bc->pStart - 1] + gfs->bc->atlasDof[gfs->bc->pEnd - gfs->bc->pStart - 1]));
1578: }
1579: gs->setup = PETSC_TRUE;
1580: PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view"));
1581: *gsection = gs;
1582: PetscFunctionReturn(PETSC_SUCCESS);
1583: }
1585: /*@
1586: PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the globallayout using
1587: a local (sequential) `PetscSection` on each MPI process and an `PetscSF` describing the section point overlap.
1589: Input Parameters:
1590: + s - The `PetscSection` for the local field layout
1591: . sf - The `PetscSF` describing parallel layout of the section points
1592: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global vector will not possess constrained dofs
1593: . numExcludes - The number of exclusion ranges, this must have the same value on all MPI processes
1594: - excludes - An array [start_0, end_0, start_1, end_1, ...] where there are `numExcludes` pairs and must have the same values on all MPI processes
1596: Output Parameter:
1597: . gsection - The `PetscSection` for the global field layout
1599: Level: advanced
1601: Notes:
1602: On each MPI process `gsection` inherits the chart of the `s` on that process.
1604: This sets negative sizes and offsets to points not owned by this process as defined by `sf` but that are within the local value of the chart of `gsection`.
1605: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1607: This routine augments `PetscSectionCreateGlobalSection()` by allowing one to exclude certain ranges in the chart of the `PetscSection`
1609: Developer Notes:
1610: This is a terrible function name
1612: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`
1613: @*/
1614: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1615: {
1616: const PetscInt *pind = NULL;
1617: PetscInt *neg = NULL, *tmpOff = NULL;
1618: PetscInt pStart, pEnd, p, e, dof, cdof, globalOff = 0, nroots;
1619: PetscInt off;
1621: PetscFunctionBegin;
1624: PetscAssertPointer(gsection, 6);
1625: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
1626: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1627: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
1628: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1629: if (nroots >= 0) {
1630: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
1631: PetscCall(PetscCalloc1(nroots, &neg));
1632: if (nroots > pEnd - pStart) {
1633: PetscCall(PetscCalloc1(nroots, &tmpOff));
1634: } else {
1635: tmpOff = &(*gsection)->atlasDof[-pStart];
1636: }
1637: }
1638: /* Mark ghost points with negative dof */
1639: for (p = pStart; p < pEnd; ++p) {
1640: for (e = 0; e < numExcludes; ++e) {
1641: if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1642: PetscCall(PetscSectionSetDof(*gsection, p, 0));
1643: break;
1644: }
1645: }
1646: if (e < numExcludes) continue;
1647: PetscCall(PetscSectionGetDof(s, p, &dof));
1648: PetscCall(PetscSectionSetDof(*gsection, p, dof));
1649: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1650: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
1651: if (neg) neg[p] = -(dof + 1);
1652: }
1653: PetscCall(PetscSectionSetUpBC(*gsection));
1654: if (nroots >= 0) {
1655: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1656: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1657: if (nroots > pEnd - pStart) {
1658: for (p = pStart; p < pEnd; ++p) {
1659: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1660: }
1661: }
1662: }
1663: /* Calculate new sizes, get process offset, and calculate point offsets */
1664: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1665: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1666: const PetscInt q = pind ? pind[p] : p;
1668: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1669: (*gsection)->atlasOff[q] = off;
1670: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1671: }
1672: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1673: globalOff -= off;
1674: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1675: (*gsection)->atlasOff[p] += globalOff;
1676: if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1677: }
1678: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1679: /* Put in negative offsets for ghost points */
1680: if (nroots >= 0) {
1681: if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1682: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1683: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1684: if (nroots > pEnd - pStart) {
1685: for (p = pStart; p < pEnd; ++p) {
1686: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1687: }
1688: }
1689: }
1690: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
1691: PetscCall(PetscFree(neg));
1692: PetscFunctionReturn(PETSC_SUCCESS);
1693: }
1695: /*@
1696: PetscSectionGetPointLayout - Get a `PetscLayout` for the points with nonzero dof counts of the unnamed default field within this `PetscSection`s local chart
1698: Collective
1700: Input Parameters:
1701: + comm - The `MPI_Comm`
1702: - s - The `PetscSection`
1704: Output Parameter:
1705: . layout - The point layout for the data that defines the section
1707: Level: advanced
1709: Notes:
1710: `PetscSectionGetValueLayout()` provides similar information but counting the total number of degrees of freedom on the MPI process (excluding constrained
1711: degrees of freedom).
1713: This count includes constrained degrees of freedom
1715: This is usually called on the default global section.
1717: Example:
1718: .vb
1719: The chart is [2,5), point 2 has 2 dof, point 3 has 0 dof, point 4 has 1 dof
1720: The local size of the `PetscLayout` is 2 since 2 points have a non-zero number of dof
1721: .ve
1723: Developer Notes:
1724: I find the names of these two functions extremely non-informative
1726: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1727: @*/
1728: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1729: {
1730: PetscInt pStart, pEnd, p, localSize = 0;
1732: PetscFunctionBegin;
1733: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1734: for (p = pStart; p < pEnd; ++p) {
1735: PetscInt dof;
1737: PetscCall(PetscSectionGetDof(s, p, &dof));
1738: if (dof >= 0) ++localSize;
1739: }
1740: PetscCall(PetscLayoutCreate(comm, layout));
1741: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1742: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1743: PetscCall(PetscLayoutSetUp(*layout));
1744: PetscFunctionReturn(PETSC_SUCCESS);
1745: }
1747: /*@
1748: PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs of a `PetscSection`
1750: Collective
1752: Input Parameters:
1753: + comm - The `MPI_Comm`
1754: - s - The `PetscSection`
1756: Output Parameter:
1757: . layout - The dof layout for the section
1759: Level: advanced
1761: Notes:
1762: `PetscSectionGetPointLayout()` provides similar information but only counting the number of points with nonzero degrees of freedom and
1763: including the constrained degrees of freedom
1765: This is usually called for the default global section.
1767: Example:
1768: .vb
1769: The chart is [2,5), point 2 has 4 dof (2 constrained), point 3 has 0 dof, point 4 has 1 dof (not constrained)
1770: The local size of the `PetscLayout` is 3 since there are 3 unconstrained degrees of freedom on this MPI process
1771: .ve
1773: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1774: @*/
1775: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1776: {
1777: PetscInt pStart, pEnd, p, localSize = 0;
1779: PetscFunctionBegin;
1781: PetscAssertPointer(layout, 3);
1782: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1783: for (p = pStart; p < pEnd; ++p) {
1784: PetscInt dof, cdof;
1786: PetscCall(PetscSectionGetDof(s, p, &dof));
1787: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1788: if (dof - cdof > 0) localSize += dof - cdof;
1789: }
1790: PetscCall(PetscLayoutCreate(comm, layout));
1791: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1792: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1793: PetscCall(PetscLayoutSetUp(*layout));
1794: PetscFunctionReturn(PETSC_SUCCESS);
1795: }
1797: /*@
1798: PetscSectionGetOffset - Return the offset into an array or `Vec` for the dof associated with the given point.
1800: Not Collective
1802: Input Parameters:
1803: + s - the `PetscSection`
1804: - point - the point
1806: Output Parameter:
1807: . offset - the offset
1809: Level: intermediate
1811: Notes:
1812: In a global section, `offset` will be negative for points not owned by this process.
1814: This is for the unnamed default field in the `PetscSection` not the named fields
1816: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1818: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetPointMajor()`
1819: @*/
1820: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1821: {
1822: PetscFunctionBegin;
1824: PetscAssertPointer(offset, 3);
1825: PetscAssert(!(point < s->pStart) && !(point >= s->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
1826: *offset = s->atlasOff[point - s->pStart];
1827: PetscFunctionReturn(PETSC_SUCCESS);
1828: }
1830: /*@
1831: PetscSectionSetOffset - Set the offset into an array or `Vec` for the dof associated with the given point.
1833: Not Collective
1835: Input Parameters:
1836: + s - the `PetscSection`
1837: . point - the point
1838: - offset - the offset, these values may be negative indicating the values are off process
1840: Level: developer
1842: Note:
1843: The user usually does not call this function, but uses `PetscSectionSetUp()`
1845: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1846: @*/
1847: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1848: {
1849: PetscFunctionBegin;
1851: PetscCheck(!(point < s->pStart) && !(point >= s->pEnd), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
1852: s->atlasOff[point - s->pStart] = offset;
1853: PetscFunctionReturn(PETSC_SUCCESS);
1854: }
1856: /*@
1857: PetscSectionGetFieldOffset - Return the offset into an array or `Vec` for the field dof associated with the given point.
1859: Not Collective
1861: Input Parameters:
1862: + s - the `PetscSection`
1863: . point - the point
1864: - field - the field
1866: Output Parameter:
1867: . offset - the offset
1869: Level: intermediate
1871: Notes:
1872: In a global section, `offset` will be negative for points not owned by this process.
1874: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1876: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldPointOffset()`
1877: @*/
1878: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1879: {
1880: PetscFunctionBegin;
1882: PetscAssertPointer(offset, 4);
1883: PetscSectionCheckValidField(field, s->numFields);
1884: PetscCall(PetscSectionGetOffset(s->field[field], point, offset));
1885: PetscFunctionReturn(PETSC_SUCCESS);
1886: }
1888: /*@
1889: PetscSectionSetFieldOffset - Set the offset into an array or `Vec` for the dof associated with the given field at a point.
1891: Not Collective
1893: Input Parameters:
1894: + s - the `PetscSection`
1895: . point - the point
1896: . field - the field
1897: - offset - the offset, these values may be negative indicating the values are off process
1899: Level: developer
1901: Note:
1902: The user usually does not call this function, but uses `PetscSectionSetUp()`
1904: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1905: @*/
1906: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1907: {
1908: PetscFunctionBegin;
1910: PetscSectionCheckValidField(field, s->numFields);
1911: PetscCall(PetscSectionSetOffset(s->field[field], point, offset));
1912: PetscFunctionReturn(PETSC_SUCCESS);
1913: }
1915: /*@
1916: PetscSectionGetFieldPointOffset - Return the offset for the first field dof associated with the given point relative to the offset for that point for the
1917: unnamed default field's first dof
1919: Not Collective
1921: Input Parameters:
1922: + s - the `PetscSection`
1923: . point - the point
1924: - field - the field
1926: Output Parameter:
1927: . offset - the offset
1929: Level: advanced
1931: Note:
1932: This ignores constraints
1934: Example:
1935: .vb
1936: if PetscSectionSetPointMajor(s,PETSC_TRUE)
1937: The unnamed default field has 3 dof at `point`
1938: Field 0 has 2 dof at `point`
1939: Then PetscSectionGetFieldPointOffset(s,point,1,&offset) returns and offset of 5
1940: .ve
1942: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldOffset()`
1943: @*/
1944: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1945: {
1946: PetscInt off, foff;
1948: PetscFunctionBegin;
1950: PetscAssertPointer(offset, 4);
1951: PetscSectionCheckValidField(field, s->numFields);
1952: PetscCall(PetscSectionGetOffset(s, point, &off));
1953: PetscCall(PetscSectionGetOffset(s->field[field], point, &foff));
1954: *offset = foff - off;
1955: PetscFunctionReturn(PETSC_SUCCESS);
1956: }
1958: /*@
1959: PetscSectionGetOffsetRange - Return the full range of offsets [`start`, `end`) for a `PetscSection`
1961: Not Collective
1963: Input Parameter:
1964: . s - the `PetscSection`
1966: Output Parameters:
1967: + start - the minimum offset
1968: - end - one more than the maximum offset
1970: Level: intermediate
1972: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1973: @*/
1974: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1975: {
1976: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1978: PetscFunctionBegin;
1980: if (s->atlasOff) {
1981: os = s->atlasOff[0];
1982: oe = s->atlasOff[0];
1983: }
1984: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1985: for (p = 0; p < pEnd - pStart; ++p) {
1986: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1988: if (off >= 0) {
1989: os = PetscMin(os, off);
1990: oe = PetscMax(oe, off + dof);
1991: }
1992: }
1993: if (start) *start = os;
1994: if (end) *end = oe;
1995: PetscFunctionReturn(PETSC_SUCCESS);
1996: }
1998: /*@
1999: PetscSectionCreateSubsection - Create a new, smaller `PetscSection` composed of only selected fields
2001: Collective
2003: Input Parameters:
2004: + s - the `PetscSection`
2005: . len - the number of subfields
2006: - fields - the subfield numbers
2008: Output Parameter:
2009: . subs - the subsection
2011: Level: advanced
2013: Notes:
2014: The chart of `subs` is the same as the chart of `s`
2016: This will error if a fieldnumber is out of range
2018: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2019: @*/
2020: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
2021: {
2022: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
2024: PetscFunctionBegin;
2025: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2027: PetscAssertPointer(fields, 3);
2028: PetscAssertPointer(subs, 4);
2029: PetscCall(PetscSectionGetNumFields(s, &nF));
2030: PetscCheck(len <= nF, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "Number of requested fields %" PetscInt_FMT " greater than number of fields %" PetscInt_FMT, len, nF);
2031: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2032: PetscCall(PetscSectionSetNumFields(*subs, len));
2033: for (f = 0; f < len; ++f) {
2034: const char *name = NULL;
2035: PetscInt numComp = 0;
2036: PetscSectionSym sym;
2038: PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
2039: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2040: PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
2041: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2042: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
2043: PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
2044: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2045: }
2046: PetscCall(PetscSectionGetFieldSym(s, fields[f], &sym));
2047: PetscCall(PetscSectionSetFieldSym(*subs, f, sym));
2048: }
2049: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2050: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2051: for (p = pStart; p < pEnd; ++p) {
2052: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
2054: for (f = 0; f < len; ++f) {
2055: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2056: PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
2057: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2058: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
2059: dof += fdof;
2060: cdof += cfdof;
2061: }
2062: PetscCall(PetscSectionSetDof(*subs, p, dof));
2063: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
2064: maxCdof = PetscMax(cdof, maxCdof);
2065: }
2066: PetscBT bst, subbst;
2068: PetscCall(PetscSectionGetBlockStarts(s, &bst));
2069: if (bst) {
2070: PetscCall(PetscBTCreate(pEnd - pStart, &subbst));
2071: PetscCall(PetscBTCopy(subbst, pEnd - pStart, bst));
2072: PetscCall(PetscSectionSetBlockStarts(*subs, subbst));
2073: }
2074: PetscCall(PetscSectionSetUp(*subs));
2075: if (maxCdof) {
2076: PetscInt *indices;
2078: PetscCall(PetscMalloc1(maxCdof, &indices));
2079: for (p = pStart; p < pEnd; ++p) {
2080: PetscInt cdof;
2082: PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
2083: if (cdof) {
2084: const PetscInt *oldIndices = NULL;
2085: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
2087: for (f = 0; f < len; ++f) {
2088: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2089: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2090: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
2091: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
2092: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
2093: numConst += cfdof;
2094: fOff += fdof;
2095: }
2096: PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
2097: }
2098: }
2099: PetscCall(PetscFree(indices));
2100: }
2101: PetscFunctionReturn(PETSC_SUCCESS);
2102: }
2104: /*@
2105: PetscSectionCreateComponentSubsection - Create a new, smaller `PetscSection` composed of only selected components
2107: Collective
2109: Input Parameters:
2110: + s - the `PetscSection`
2111: . len - the number of components
2112: - comps - the component numbers
2114: Output Parameter:
2115: . subs - the subsection
2117: Level: advanced
2119: Notes:
2120: The chart of `subs` is the same as the chart of `s`
2122: This will error if the section has more than one field, or if a component number is out of range
2124: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2125: @*/
2126: PetscErrorCode PetscSectionCreateComponentSubsection(PetscSection s, PetscInt len, const PetscInt comps[], PetscSection *subs)
2127: {
2128: PetscSectionSym sym;
2129: const char *name = NULL;
2130: PetscInt Nf, pStart, pEnd;
2132: PetscFunctionBegin;
2133: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2135: PetscAssertPointer(comps, 3);
2136: PetscAssertPointer(subs, 4);
2137: PetscCall(PetscSectionGetNumFields(s, &Nf));
2138: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "This method can only handle one field, not %" PetscInt_FMT, Nf);
2139: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2140: PetscCall(PetscSectionSetNumFields(*subs, 1));
2141: PetscCall(PetscSectionGetFieldName(s, 0, &name));
2142: PetscCall(PetscSectionSetFieldName(*subs, 0, name));
2143: PetscCall(PetscSectionSetFieldComponents(*subs, 0, len));
2144: PetscCall(PetscSectionGetFieldSym(s, 0, &sym));
2145: PetscCall(PetscSectionSetFieldSym(*subs, 0, sym));
2146: for (PetscInt c = 0; c < len; ++c) {
2147: PetscCall(PetscSectionGetComponentName(s, 0, comps[c], &name));
2148: PetscCall(PetscSectionSetComponentName(*subs, 0, c, name));
2149: }
2150: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2151: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2152: for (PetscInt p = pStart; p < pEnd; ++p) {
2153: PetscInt dof, cdof, cfdof;
2155: PetscCall(PetscSectionGetDof(s, p, &dof));
2156: if (!dof) continue;
2157: PetscCall(PetscSectionGetFieldConstraintDof(s, p, 0, &cfdof));
2158: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2159: PetscCheck(!cdof && !cfdof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Component selection does not work with constraints");
2160: PetscCall(PetscSectionSetFieldDof(*subs, p, 0, len));
2161: PetscCall(PetscSectionSetDof(*subs, p, len));
2162: }
2163: PetscCall(PetscSectionSetUp(*subs));
2164: PetscFunctionReturn(PETSC_SUCCESS);
2165: }
2167: /*@
2168: PetscSectionCreateSupersection - Create a new, larger section composed of multiple `PetscSection`s
2170: Collective
2172: Input Parameters:
2173: + s - the input sections
2174: - len - the number of input sections
2176: Output Parameter:
2177: . supers - the supersection
2179: Level: advanced
2181: Notes:
2182: The section offsets now refer to a new, larger vector.
2184: Developer Notes:
2185: Needs to explain how the sections are composed
2187: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2188: @*/
2189: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2190: {
2191: PetscInt Nf = 0, f, pStart = PETSC_INT_MAX, pEnd = 0, p, maxCdof = 0, i;
2193: PetscFunctionBegin;
2194: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2195: for (i = 0; i < len; ++i) {
2196: PetscInt nf, pStarti, pEndi;
2198: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2199: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2200: pStart = PetscMin(pStart, pStarti);
2201: pEnd = PetscMax(pEnd, pEndi);
2202: Nf += nf;
2203: }
2204: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2205: PetscCall(PetscSectionSetNumFields(*supers, Nf));
2206: for (i = 0, f = 0; i < len; ++i) {
2207: PetscInt nf, fi, ci;
2209: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2210: for (fi = 0; fi < nf; ++fi, ++f) {
2211: const char *name = NULL;
2212: PetscInt numComp = 0;
2214: PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2215: PetscCall(PetscSectionSetFieldName(*supers, f, name));
2216: PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2217: PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2218: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2219: PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2220: PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2221: }
2222: }
2223: }
2224: PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2225: for (p = pStart; p < pEnd; ++p) {
2226: PetscInt dof = 0, cdof = 0;
2228: for (i = 0, f = 0; i < len; ++i) {
2229: PetscInt nf, fi, pStarti, pEndi;
2230: PetscInt fdof = 0, cfdof = 0;
2232: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2233: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2234: if ((p < pStarti) || (p >= pEndi)) continue;
2235: for (fi = 0; fi < nf; ++fi, ++f) {
2236: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2237: PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2238: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2239: if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2240: dof += fdof;
2241: cdof += cfdof;
2242: }
2243: }
2244: PetscCall(PetscSectionSetDof(*supers, p, dof));
2245: if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2246: maxCdof = PetscMax(cdof, maxCdof);
2247: }
2248: PetscCall(PetscSectionSetUp(*supers));
2249: if (maxCdof) {
2250: PetscInt *indices;
2252: PetscCall(PetscMalloc1(maxCdof, &indices));
2253: for (p = pStart; p < pEnd; ++p) {
2254: PetscInt cdof;
2256: PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2257: if (cdof) {
2258: PetscInt dof, numConst = 0, fOff = 0;
2260: for (i = 0, f = 0; i < len; ++i) {
2261: const PetscInt *oldIndices = NULL;
2262: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
2264: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2265: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2266: if ((p < pStarti) || (p >= pEndi)) continue;
2267: for (fi = 0; fi < nf; ++fi, ++f) {
2268: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2269: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2270: PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2271: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2272: PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2273: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2274: numConst += cfdof;
2275: }
2276: PetscCall(PetscSectionGetDof(s[i], p, &dof));
2277: fOff += dof;
2278: }
2279: PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2280: }
2281: }
2282: PetscCall(PetscFree(indices));
2283: }
2284: PetscFunctionReturn(PETSC_SUCCESS);
2285: }
2287: static PetscErrorCode PetscSectionCreateSubplexSection_Private(PetscSection s, IS subpointIS, PetscBool renumberPoints, PetscSection *subs)
2288: {
2289: const PetscInt *points = NULL, *indices = NULL;
2290: PetscInt *spoints = NULL, *order = NULL;
2291: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
2293: PetscFunctionBegin;
2296: PetscAssertPointer(subs, 4);
2297: PetscCall(PetscSectionGetNumFields(s, &numFields));
2298: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2299: if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2300: for (f = 0; f < numFields; ++f) {
2301: const char *name = NULL;
2302: PetscInt numComp = 0;
2304: PetscCall(PetscSectionGetFieldName(s, f, &name));
2305: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2306: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2307: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2308: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2309: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2310: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2311: }
2312: }
2313: /* For right now, we do not try to squeeze the subchart */
2314: if (subpointIS) {
2315: PetscCall(ISGetLocalSize(subpointIS, &numSubpoints));
2316: PetscCall(ISGetIndices(subpointIS, &points));
2317: }
2318: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2319: if (renumberPoints) {
2320: PetscBool sorted;
2322: spStart = 0;
2323: spEnd = numSubpoints;
2324: PetscCall(ISSorted(subpointIS, &sorted));
2325: if (!sorted) {
2326: PetscCall(PetscMalloc2(numSubpoints, &spoints, numSubpoints, &order));
2327: PetscCall(PetscArraycpy(spoints, points, numSubpoints));
2328: for (PetscInt i = 0; i < numSubpoints; ++i) order[i] = i;
2329: PetscCall(PetscSortIntWithArray(numSubpoints, spoints, order));
2330: }
2331: } else {
2332: PetscCall(ISGetMinMax(subpointIS, &spStart, &spEnd));
2333: ++spEnd;
2334: }
2335: PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2336: for (p = pStart; p < pEnd; ++p) {
2337: PetscInt dof, cdof, fdof = 0, cfdof = 0;
2339: PetscCall(PetscFindInt(p, numSubpoints, spoints ? spoints : points, &subp));
2340: if (subp < 0) continue;
2341: if (!renumberPoints) subp = p;
2342: else subp = order ? order[subp] : subp;
2343: for (f = 0; f < numFields; ++f) {
2344: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2345: PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2346: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2347: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2348: }
2349: PetscCall(PetscSectionGetDof(s, p, &dof));
2350: PetscCall(PetscSectionSetDof(*subs, subp, dof));
2351: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2352: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2353: }
2354: PetscCall(PetscSectionSetUp(*subs));
2355: /* Change offsets to original offsets */
2356: for (p = pStart; p < pEnd; ++p) {
2357: PetscInt off, foff = 0;
2359: PetscCall(PetscFindInt(p, numSubpoints, spoints ? spoints : points, &subp));
2360: if (subp < 0) continue;
2361: if (!renumberPoints) subp = p;
2362: else subp = order ? order[subp] : subp;
2363: for (f = 0; f < numFields; ++f) {
2364: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2365: PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2366: }
2367: PetscCall(PetscSectionGetOffset(s, p, &off));
2368: PetscCall(PetscSectionSetOffset(*subs, subp, off));
2369: }
2370: /* Copy constraint indices */
2371: for (subp = spStart; subp < spEnd; ++subp) {
2372: PetscInt cdof;
2374: PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2375: if (cdof) {
2376: for (f = 0; f < numFields; ++f) {
2377: PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2378: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2379: }
2380: PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2381: PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2382: }
2383: }
2384: if (subpointIS) PetscCall(ISRestoreIndices(subpointIS, &points));
2385: PetscCall(PetscFree2(spoints, order));
2386: PetscFunctionReturn(PETSC_SUCCESS);
2387: }
2389: /*@
2390: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2392: Collective
2394: Input Parameters:
2395: + s - the `PetscSection`
2396: - subpointIS - a sorted list of points in the original mesh which are in the submesh
2398: Output Parameter:
2399: . subs - the subsection
2401: Level: advanced
2403: Notes:
2404: The points are renumbered from 0, and the section offsets now refer to a new, smaller vector. That is the chart of `subs` is `[0,sizeof(subpointmap))`
2406: Compare this with `PetscSectionCreateSubdomainSection()` that does not map the points numbers to start at zero but leaves them as before
2408: Developer Notes:
2409: The use of the term Submesh is confusing and needs clarification, it is not specific to meshes. It appears to be just a subset of the chart of the original `PetscSection`
2411: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2412: @*/
2413: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointIS, PetscSection *subs)
2414: {
2415: PetscFunctionBegin;
2416: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointIS, PETSC_TRUE, subs));
2417: PetscFunctionReturn(PETSC_SUCCESS);
2418: }
2420: /*@
2421: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2423: Collective
2425: Input Parameters:
2426: + s - the `PetscSection`
2427: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2429: Output Parameter:
2430: . subs - the subsection
2432: Level: advanced
2434: Notes:
2435: The point numbers remain the same as in the larger `PetscSection`, but the section offsets now refer to a new, smaller vector. The chart of `subs`
2436: is `[min(subpointMap),max(subpointMap)+1)`
2438: Compare this with `PetscSectionCreateSubmeshSection()` that maps the point numbers to start at zero
2440: Developer Notes:
2441: The use of the term Subdomain is unneeded and needs clarification, it is not specific to meshes. It appears to be just a subset of the chart of the original `PetscSection`
2443: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2444: @*/
2445: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2446: {
2447: PetscFunctionBegin;
2448: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_FALSE, subs));
2449: PetscFunctionReturn(PETSC_SUCCESS);
2450: }
2452: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2453: {
2454: PetscMPIInt rank;
2456: PetscFunctionBegin;
2457: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2458: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2459: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2460: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
2461: if (s->bc && s->bc->atlasDof[p] > 0) {
2462: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2463: if (s->bcIndices) {
2464: for (PetscInt b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2465: }
2466: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2467: } else {
2468: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2469: }
2470: }
2471: PetscCall(PetscViewerFlush(viewer));
2472: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2473: if (s->sym) {
2474: PetscCall(PetscViewerASCIIPushTab(viewer));
2475: PetscCall(PetscSectionSymView(s->sym, viewer));
2476: PetscCall(PetscViewerASCIIPopTab(viewer));
2477: }
2478: PetscFunctionReturn(PETSC_SUCCESS);
2479: }
2481: /*@
2482: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2484: Collective
2486: Input Parameters:
2487: + A - the `PetscSection` object to view
2488: . obj - optional object that provides the options prefix used for the options, pass `NULL` to use the options prefix of `A`
2489: - name - command line option
2491: Options Database Key:
2492: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
2494: Level: intermediate
2496: Note:
2497: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
2498: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
2500: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionView()`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscOptionsCreateViewer()`
2501: @*/
2502: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2503: {
2504: PetscFunctionBegin;
2506: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2507: PetscFunctionReturn(PETSC_SUCCESS);
2508: }
2510: /*@
2511: PetscSectionView - Views a `PetscSection`
2513: Collective
2515: Input Parameters:
2516: + s - the `PetscSection` object to view
2517: - viewer - the viewer
2519: Level: beginner
2521: Note:
2522: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2523: distribution independent data, such as dofs, offsets, constraint dofs,
2524: and constraint indices. Points that have negative dofs, for instance,
2525: are not saved as they represent points owned by other processes.
2526: Point numbering and rank assignment is currently not stored.
2527: The saved section can be loaded with `PetscSectionLoad()`.
2529: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2530: @*/
2531: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2532: {
2533: PetscBool isascii, ishdf5;
2534: PetscInt f;
2536: PetscFunctionBegin;
2538: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2540: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2541: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2542: if (isascii) {
2543: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2544: if (s->numFields) {
2545: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2546: for (f = 0; f < s->numFields; ++f) {
2547: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " \"%s\" with %" PetscInt_FMT " components\n", f, s->fieldNames[f], s->numFieldComponents[f]));
2548: PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2549: }
2550: } else {
2551: PetscCall(PetscSectionView_ASCII(s, viewer));
2552: }
2553: } else if (ishdf5) {
2554: #if PetscDefined(HAVE_HDF5)
2555: PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2556: #else
2557: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2558: #endif
2559: }
2560: PetscFunctionReturn(PETSC_SUCCESS);
2561: }
2563: /*@
2564: PetscSectionLoad - Loads a `PetscSection`
2566: Collective
2568: Input Parameters:
2569: + s - the `PetscSection` object to load
2570: - viewer - the viewer
2572: Level: beginner
2574: Note:
2575: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2576: a section saved with `PetscSectionView()`. The number of processes
2577: used here (N) does not need to be the same as that used when saving.
2578: After calling this function, the chart of s on rank i will be set
2579: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2580: saved section points.
2582: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2583: @*/
2584: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2585: {
2586: PetscBool ishdf5;
2588: PetscFunctionBegin;
2591: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2592: PetscCheck(ishdf5, PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2593: #if PetscDefined(HAVE_HDF5)
2594: PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2595: PetscFunctionReturn(PETSC_SUCCESS);
2596: #else
2597: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2598: #endif
2599: }
2601: static inline PetscErrorCode PrintArrayElement(void *array, PetscDataType data_type, PetscCount index, PetscViewer viewer)
2602: {
2603: PetscFunctionBeginUser;
2604: switch (data_type) {
2605: case PETSC_INT: {
2606: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt_FMT, ((PetscInt *)array)[index]));
2607: break;
2608: }
2609: case PETSC_INT32: {
2610: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt32_FMT, ((PetscInt32 *)array)[index]));
2611: break;
2612: }
2613: case PETSC_INT64: {
2614: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt64_FMT, ((PetscInt64 *)array)[index]));
2615: break;
2616: }
2617: case PETSC_COUNT: {
2618: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscCount_FMT, ((PetscCount *)array)[index]));
2619: break;
2620: }
2621: // PETSC_SCALAR is set to the appropriate type
2622: case PETSC_DOUBLE: {
2623: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", ((double *)array)[index]));
2624: break;
2625: }
2626: case PETSC_FLOAT: {
2627: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((float *)array)[index]));
2628: break;
2629: }
2630: #if PetscDefined(USE_REAL___FLOAT128)
2631: case PETSC___FLOAT128: {
2632: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((PetscReal *)array)[index]));
2633: break;
2634: }
2635: #endif
2636: #if PetscDefined(USE_REAL___FP16)
2637: case PETSC___FP16: {
2638: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((PetscReal *)array)[index]));
2639: break;
2640: }
2641: #endif
2642: #if PetscDefined(HAVE_COMPLEX)
2643: case PETSC_COMPLEX: {
2644: PetscComplex v = ((PetscComplex *)array)[index];
2645: if (PetscImaginaryPartComplex(v) > 0.0) {
2646: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPartComplex(v), (double)PetscImaginaryPartComplex(v)));
2647: } else if (PetscImaginaryPartComplex(v) < 0.0) {
2648: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPartComplex(v), (double)(-PetscImaginaryPartComplex(v))));
2649: } else {
2650: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPartComplex(v)));
2651: }
2652: break;
2653: }
2654: #endif
2655: default:
2656: SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "PetscDataType %d (%s) not supported", data_type, PetscDataTypes[data_type]);
2657: }
2658: PetscFunctionReturn(PETSC_SUCCESS);
2659: }
2661: PetscErrorCode PetscSectionArrayView_ASCII_Internal(PetscSection s, void *array, PetscDataType data_type, PetscViewer viewer)
2662: {
2663: PetscInt i;
2664: PetscMPIInt rank;
2666: PetscFunctionBegin;
2667: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2668: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2669: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2670: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
2671: if (s->bc && (s->bc->atlasDof[p] > 0)) {
2672: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2673: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) PetscCall(PrintArrayElement(array, data_type, i, viewer));
2674: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " constrained"));
2675: for (PetscInt b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2676: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2677: } else {
2678: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2679: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) PetscCall(PrintArrayElement(array, data_type, i, viewer));
2680: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2681: }
2682: }
2683: PetscCall(PetscViewerFlush(viewer));
2684: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2685: PetscFunctionReturn(PETSC_SUCCESS);
2686: }
2688: /*@
2689: PetscSectionArrayView - View an array, using the section to structure the values
2691: Collective
2693: Input Parameters:
2694: + s - the organizing `PetscSection`
2695: . array - the array of values
2696: . data_type - the `PetscDataType` of the array
2697: - viewer - the `PetscViewer`
2699: Level: developer
2701: .seealso: `PetscSection`, `PetscViewer`, `PetscSectionCreate()`, `VecSetValuesSection()`, `PetscSectionVecView()`
2702: @*/
2703: PetscErrorCode PetscSectionArrayView(PetscSection s, void *array, PetscDataType data_type, PetscViewer viewer)
2704: {
2705: PetscBool isascii;
2706: PetscInt f;
2708: PetscFunctionBegin;
2710: if (!array) {
2711: PetscInt size;
2712: PetscCall(PetscSectionGetStorageSize(s, &size));
2713: PetscCheck(size == 0, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_SIZ, "NULL array passed, but section's storage size is non-zero");
2714: } else PetscAssertPointer(array, 2);
2715: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2717: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2718: if (isascii) {
2719: if (s->numFields) {
2720: PetscCall(PetscViewerASCIIPrintf(viewer, "Array with %" PetscInt_FMT " fields\n", s->numFields));
2721: for (f = 0; f < s->numFields; ++f) {
2722: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
2723: PetscCall(PetscSectionArrayView_ASCII_Internal(s->field[f], array, data_type, viewer));
2724: }
2725: } else {
2726: PetscCall(PetscSectionArrayView_ASCII_Internal(s, array, data_type, viewer));
2727: }
2728: }
2729: PetscFunctionReturn(PETSC_SUCCESS);
2730: }
2732: /*@
2733: PetscSectionResetClosurePermutation - Remove any existing closure permutation
2735: Input Parameter:
2736: . section - The `PetscSection`
2738: Level: intermediate
2740: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionSetClosureIndex()`, `PetscSectionReset()`
2741: @*/
2742: PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2743: {
2744: PetscSectionClosurePermVal clVal;
2746: PetscFunctionBegin;
2747: if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2748: kh_foreach_value(section->clHash, clVal, {
2749: PetscCall(PetscFree(clVal.perm));
2750: PetscCall(PetscFree(clVal.invPerm));
2751: });
2752: kh_destroy(ClPerm, section->clHash);
2753: section->clHash = NULL;
2754: PetscFunctionReturn(PETSC_SUCCESS);
2755: }
2757: /*@
2758: PetscSectionReset - Frees all section data, the section is then as if `PetscSectionCreate()` had just been called.
2760: Not Collective
2762: Input Parameter:
2763: . s - the `PetscSection`
2765: Level: beginner
2767: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`
2768: @*/
2769: PetscErrorCode PetscSectionReset(PetscSection s)
2770: {
2771: PetscInt f, c;
2773: PetscFunctionBegin;
2775: for (f = 0; f < s->numFields; ++f) {
2776: PetscCall(PetscSectionDestroy(&s->field[f]));
2777: PetscCall(PetscFree(s->fieldNames[f]));
2778: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2779: PetscCall(PetscFree(s->compNames[f]));
2780: }
2781: PetscCall(PetscFree(s->numFieldComponents));
2782: PetscCall(PetscFree(s->fieldNames));
2783: PetscCall(PetscFree(s->compNames));
2784: PetscCall(PetscFree(s->field));
2785: PetscCall(PetscSectionDestroy(&s->bc));
2786: PetscCall(PetscFree(s->bcIndices));
2787: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2788: PetscCall(PetscSectionDestroy(&s->clSection));
2789: PetscCall(ISDestroy(&s->clPoints));
2790: PetscCall(ISDestroy(&s->perm));
2791: PetscCall(PetscBTDestroy(&s->blockStarts));
2792: PetscCall(PetscSectionResetClosurePermutation(s));
2793: PetscCall(PetscSectionSymDestroy(&s->sym));
2794: PetscCall(PetscSectionDestroy(&s->clSection));
2795: PetscCall(ISDestroy(&s->clPoints));
2796: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2797: s->pStart = -1;
2798: s->pEnd = -1;
2799: s->maxDof = 0;
2800: s->setup = PETSC_FALSE;
2801: s->numFields = 0;
2802: s->clObj = NULL;
2803: PetscFunctionReturn(PETSC_SUCCESS);
2804: }
2806: /*@
2807: PetscSectionDestroy - Frees a `PetscSection`
2809: Not Collective
2811: Input Parameter:
2812: . s - the `PetscSection`
2814: Level: beginner
2816: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2817: @*/
2818: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2819: {
2820: PetscFunctionBegin;
2821: if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2823: if (--((PetscObject)*s)->refct > 0) {
2824: *s = NULL;
2825: PetscFunctionReturn(PETSC_SUCCESS);
2826: }
2827: PetscCall(PetscSectionReset(*s));
2828: PetscCall(PetscHeaderDestroy(s));
2829: PetscFunctionReturn(PETSC_SUCCESS);
2830: }
2832: static PetscErrorCode VecIntGetValuesSection_Private(const PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2833: {
2834: const PetscInt p = point - s->pStart;
2836: PetscFunctionBegin;
2838: *values = &baseArray[s->atlasOff[p]];
2839: PetscFunctionReturn(PETSC_SUCCESS);
2840: }
2842: static PetscErrorCode VecIntSetValuesSection_Private(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2843: {
2844: PetscInt *array;
2845: const PetscInt p = point - s->pStart;
2846: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2847: PetscInt cDim = 0;
2849: PetscFunctionBegin;
2851: PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2852: array = &baseArray[s->atlasOff[p]];
2853: if (!cDim) {
2854: if (orientation >= 0) {
2855: const PetscInt dim = s->atlasDof[p];
2856: PetscInt i;
2858: if (mode == INSERT_VALUES) {
2859: for (i = 0; i < dim; ++i) array[i] = values ? values[i] : i;
2860: } else {
2861: for (i = 0; i < dim; ++i) array[i] += values[i];
2862: }
2863: } else {
2864: PetscInt offset = 0;
2865: PetscInt j = -1, field, i;
2867: for (field = 0; field < s->numFields; ++field) {
2868: const PetscInt dim = s->field[field]->atlasDof[p];
2870: for (i = dim - 1; i >= 0; --i) array[++j] = values ? values[i + offset] : i + offset;
2871: offset += dim;
2872: }
2873: }
2874: } else {
2875: if (orientation >= 0) {
2876: const PetscInt dim = s->atlasDof[p];
2877: PetscInt cInd = 0, i;
2878: const PetscInt *cDof;
2880: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2881: if (mode == INSERT_VALUES) {
2882: for (i = 0; i < dim; ++i) {
2883: if ((cInd < cDim) && (i == cDof[cInd])) {
2884: ++cInd;
2885: continue;
2886: }
2887: array[i] = values ? values[i] : i;
2888: }
2889: } else {
2890: for (i = 0; i < dim; ++i) {
2891: if ((cInd < cDim) && (i == cDof[cInd])) {
2892: ++cInd;
2893: continue;
2894: }
2895: array[i] += values[i];
2896: }
2897: }
2898: } else {
2899: const PetscInt *cDof;
2900: PetscInt offset = 0;
2901: PetscInt cOffset = 0;
2902: PetscInt j = 0, field;
2904: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2905: for (field = 0; field < s->numFields; ++field) {
2906: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2907: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2908: const PetscInt sDim = dim - tDim;
2909: PetscInt cInd = 0, i, k;
2911: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2912: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2913: ++cInd;
2914: continue;
2915: }
2916: array[j] = values ? values[k] : k;
2917: }
2918: offset += dim;
2919: cOffset += dim - tDim;
2920: }
2921: }
2922: }
2923: PetscFunctionReturn(PETSC_SUCCESS);
2924: }
2926: /*@
2927: PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs
2929: Not Collective
2931: Input Parameter:
2932: . s - The `PetscSection`
2934: Output Parameter:
2935: . hasConstraints - flag indicating that the section has constrained dofs
2937: Level: intermediate
2939: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2940: @*/
2941: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2942: {
2943: PetscFunctionBegin;
2945: PetscAssertPointer(hasConstraints, 2);
2946: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2947: PetscFunctionReturn(PETSC_SUCCESS);
2948: }
2950: /*@
2951: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained for a given point
2953: Not Collective
2955: Input Parameters:
2956: + s - The `PetscSection`
2957: - point - The point
2959: Output Parameter:
2960: . indices - The constrained dofs
2962: Level: intermediate
2964: Fortran Notes:
2965: Use `PetscSectionRestoreConstraintIndices()` when the indices are no longer needed
2967: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2968: @*/
2969: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt *indices[])
2970: {
2971: PetscFunctionBegin;
2973: if (s->bc) PetscCall(VecIntGetValuesSection_Private(s->bcIndices, s->bc, point, indices));
2974: else *indices = NULL;
2975: PetscFunctionReturn(PETSC_SUCCESS);
2976: }
2978: /*@
2979: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2981: Not Collective
2983: Input Parameters:
2984: + s - The `PetscSection`
2985: . point - The point
2986: - indices - The constrained dofs
2988: Level: intermediate
2990: .seealso: [PetscSection](ch_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2991: @*/
2992: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2993: {
2994: PetscFunctionBegin;
2996: if (s->bc) {
2997: const PetscInt dof = s->atlasDof[point];
2998: const PetscInt cdof = s->bc->atlasDof[point];
2999: if (indices)
3000: for (PetscInt d = 0; d < cdof; ++d)
3001: PetscCheck(indices[d] < dof, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " dof %" PetscInt_FMT ", invalid constraint index[%" PetscInt_FMT "]: %" PetscInt_FMT, point, dof, d, indices[d]);
3002: PetscCall(VecIntSetValuesSection_Private(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
3003: }
3004: PetscFunctionReturn(PETSC_SUCCESS);
3005: }
3007: /*@
3008: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
3010: Not Collective
3012: Input Parameters:
3013: + s - The `PetscSection`
3014: . field - The field number
3015: - point - The point
3017: Output Parameter:
3018: . indices - The constrained dofs sorted in ascending order, the length is returned by `PetscSectionGetConstraintDof()`.
3020: Level: intermediate
3022: Fortran Notes:
3023: Use `PetscSectionRestoreFieldConstraintIndices()` to restore the indices when no longer needed
3025: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
3026: @*/
3027: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt *indices[])
3028: {
3029: PetscFunctionBegin;
3031: PetscAssertPointer(indices, 4);
3032: PetscSectionCheckValidField(field, s->numFields);
3033: PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
3034: PetscFunctionReturn(PETSC_SUCCESS);
3035: }
3037: /*@
3038: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
3040: Not Collective
3042: Input Parameters:
3043: + s - The `PetscSection`
3044: . point - The point
3045: . field - The field number
3046: - indices - The constrained dofs
3048: Level: intermediate
3050: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
3051: @*/
3052: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
3053: {
3054: PetscFunctionBegin;
3056: PetscSectionCheckValidField(field, s->numFields);
3057: PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
3058: PetscFunctionReturn(PETSC_SUCCESS);
3059: }
3061: /*@
3062: PetscSectionPermute - Reorder the section according to the input point permutation
3064: Collective
3066: Input Parameters:
3067: + section - The `PetscSection` object
3068: - permutation - The point permutation, old point p becomes new point perm[p]
3070: Output Parameter:
3071: . sectionNew - The permuted `PetscSection`
3073: Level: intermediate
3075: Note:
3076: The data and the access to the data via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` are both changed in `sectionNew`
3078: Compare to `PetscSectionSetPermutation()`
3080: .seealso: [PetscSection](ch_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
3081: @*/
3082: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
3083: {
3084: PetscSection s = section, sNew;
3085: const PetscInt *perm;
3086: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
3088: PetscFunctionBegin;
3091: PetscAssertPointer(sectionNew, 3);
3092: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
3093: PetscCall(PetscSectionGetNumFields(s, &numFields));
3094: if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
3095: for (f = 0; f < numFields; ++f) {
3096: const char *name;
3097: PetscInt numComp;
3099: PetscCall(PetscSectionGetFieldName(s, f, &name));
3100: PetscCall(PetscSectionSetFieldName(sNew, f, name));
3101: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
3102: PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
3103: for (c = 0; c < s->numFieldComponents[f]; ++c) {
3104: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
3105: PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
3106: }
3107: }
3108: PetscCall(ISGetLocalSize(permutation, &numPoints));
3109: PetscCall(ISGetIndices(permutation, &perm));
3110: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
3111: PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
3112: PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
3113: for (p = pStart; p < pEnd; ++p) {
3114: PetscInt dof, cdof;
3116: PetscCall(PetscSectionGetDof(s, p, &dof));
3117: PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
3118: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
3119: if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
3120: for (f = 0; f < numFields; ++f) {
3121: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
3122: PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
3123: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
3124: if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
3125: }
3126: }
3127: PetscCall(PetscSectionSetUp(sNew));
3128: for (p = pStart; p < pEnd; ++p) {
3129: const PetscInt *cind;
3130: PetscInt cdof;
3132: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
3133: if (cdof) {
3134: PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
3135: PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
3136: }
3137: for (f = 0; f < numFields; ++f) {
3138: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
3139: if (cdof) {
3140: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
3141: PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
3142: }
3143: }
3144: }
3145: PetscCall(ISRestoreIndices(permutation, &perm));
3146: *sectionNew = sNew;
3147: PetscFunctionReturn(PETSC_SUCCESS);
3148: }
3150: /*@
3151: PetscSectionSetClosureIndex - Create an internal data structure to speed up closure queries.
3153: Collective
3155: Input Parameters:
3156: + section - The `PetscSection`
3157: . obj - A `PetscObject` which serves as the key for this index
3158: . clSection - `PetscSection` giving the size of the closure of each point
3159: - clPoints - `IS` giving the points in each closure
3161: Level: advanced
3163: Note:
3164: This function creates an internal map from each point to its closure. We compress out closure points with no dofs in this section.
3166: Developer Notes:
3167: The information provided here is completely opaque
3169: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
3170: @*/
3171: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
3172: {
3173: PetscFunctionBegin;
3177: if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
3178: section->clObj = obj;
3179: PetscCall(PetscObjectReference((PetscObject)clSection));
3180: PetscCall(PetscObjectReference((PetscObject)clPoints));
3181: PetscCall(PetscSectionDestroy(§ion->clSection));
3182: PetscCall(ISDestroy(§ion->clPoints));
3183: section->clSection = clSection;
3184: section->clPoints = clPoints;
3185: PetscFunctionReturn(PETSC_SUCCESS);
3186: }
3188: /*@
3189: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
3191: Collective
3193: Input Parameters:
3194: + section - The `PetscSection`
3195: - obj - A `PetscObject` which serves as the key for this index
3197: Output Parameters:
3198: + clSection - `PetscSection` giving the size of the closure of each point
3199: - clPoints - `IS` giving the points in each closure
3201: Level: advanced
3203: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3204: @*/
3205: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
3206: {
3207: PetscFunctionBegin;
3208: if (section->clObj == obj) {
3209: if (clSection) *clSection = section->clSection;
3210: if (clPoints) *clPoints = section->clPoints;
3211: } else {
3212: if (clSection) *clSection = NULL;
3213: if (clPoints) *clPoints = NULL;
3214: }
3215: PetscFunctionReturn(PETSC_SUCCESS);
3216: }
3218: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
3219: {
3220: khiter_t iter;
3221: int new_entry;
3222: PetscSectionClosurePermKey key = {depth, clSize};
3223: PetscSectionClosurePermVal *val;
3225: PetscFunctionBegin;
3226: if (section->clObj != obj) {
3227: PetscCall(PetscSectionDestroy(§ion->clSection));
3228: PetscCall(ISDestroy(§ion->clPoints));
3229: }
3230: section->clObj = obj;
3231: if (!section->clHash) PetscCall(PetscClPermCreate(§ion->clHash));
3232: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
3233: val = &kh_val(section->clHash, iter);
3234: if (!new_entry) {
3235: PetscCall(PetscFree(val->perm));
3236: PetscCall(PetscFree(val->invPerm));
3237: }
3238: if (mode == PETSC_COPY_VALUES) {
3239: PetscCall(PetscMalloc1(clSize, &val->perm));
3240: PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
3241: } else if (mode == PETSC_OWN_POINTER) {
3242: val->perm = clPerm;
3243: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
3244: PetscCall(PetscMalloc1(clSize, &val->invPerm));
3245: for (PetscInt i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
3246: PetscFunctionReturn(PETSC_SUCCESS);
3247: }
3249: /*@
3250: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3252: Not Collective
3254: Input Parameters:
3255: + section - The `PetscSection`
3256: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3257: . depth - Depth of points on which to apply the given permutation
3258: - perm - Permutation of the cell dof closure
3260: Level: intermediate
3262: Notes:
3263: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
3264: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
3265: each topology and degree.
3267: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
3268: exotic/enriched spaces on mixed topology meshes.
3270: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
3271: @*/
3272: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
3273: {
3274: const PetscInt *clPerm = NULL;
3275: PetscInt clSize = 0;
3277: PetscFunctionBegin;
3278: if (perm) {
3279: PetscCall(ISGetLocalSize(perm, &clSize));
3280: PetscCall(ISGetIndices(perm, &clPerm));
3281: }
3282: PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
3283: if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
3284: PetscFunctionReturn(PETSC_SUCCESS);
3285: }
3287: static PetscErrorCode PetscSectionGetClosurePermutation_Private(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3288: {
3289: PetscFunctionBegin;
3290: if (section->clObj == obj) {
3291: PetscSectionClosurePermKey k = {depth, size};
3292: PetscSectionClosurePermVal v;
3294: PetscCall(PetscClPermGet(section->clHash, k, &v));
3295: if (perm) *perm = v.perm;
3296: } else {
3297: if (perm) *perm = NULL;
3298: }
3299: PetscFunctionReturn(PETSC_SUCCESS);
3300: }
3302: /*@
3303: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3305: Not Collective
3307: Input Parameters:
3308: + section - The `PetscSection`
3309: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
3310: . depth - Depth stratum on which to obtain closure permutation
3311: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3313: Output Parameter:
3314: . perm - The dof closure permutation
3316: Level: intermediate
3318: Note:
3319: The user must destroy the `IS` that is returned.
3321: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3322: @*/
3323: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3324: {
3325: const PetscInt *clPerm = NULL;
3327: PetscFunctionBegin;
3328: PetscCall(PetscSectionGetClosurePermutation_Private(section, obj, depth, clSize, &clPerm));
3329: PetscCheck(clPerm, PetscObjectComm(obj), PETSC_ERR_ARG_WRONG, "There is no closure permutation associated with this object for depth %" PetscInt_FMT " of size %" PetscInt_FMT, depth, clSize);
3330: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3331: PetscFunctionReturn(PETSC_SUCCESS);
3332: }
3334: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3335: {
3336: PetscFunctionBegin;
3337: if (section->clObj == obj && section->clHash) {
3338: PetscSectionClosurePermKey k = {depth, size};
3339: PetscSectionClosurePermVal v;
3340: PetscCall(PetscClPermGet(section->clHash, k, &v));
3341: if (perm) *perm = v.invPerm;
3342: } else {
3343: if (perm) *perm = NULL;
3344: }
3345: PetscFunctionReturn(PETSC_SUCCESS);
3346: }
3348: /*@
3349: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
3351: Not Collective
3353: Input Parameters:
3354: + section - The `PetscSection`
3355: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3356: . depth - Depth stratum on which to obtain closure permutation
3357: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3359: Output Parameter:
3360: . perm - The dof closure permutation
3362: Level: intermediate
3364: Note:
3365: The user must destroy the `IS` that is returned.
3367: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3368: @*/
3369: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3370: {
3371: const PetscInt *clPerm = NULL;
3373: PetscFunctionBegin;
3374: PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3375: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3376: PetscFunctionReturn(PETSC_SUCCESS);
3377: }
3379: /*@
3380: PetscSectionGetField - Get the `PetscSection` associated with a single field
3382: Input Parameters:
3383: + s - The `PetscSection`
3384: - field - The field number
3386: Output Parameter:
3387: . subs - The `PetscSection` for the given field, note the chart of `subs` is not set
3389: Level: intermediate
3391: Note:
3392: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
3394: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3395: @*/
3396: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3397: {
3398: PetscFunctionBegin;
3400: PetscAssertPointer(subs, 3);
3401: PetscSectionCheckValidField(field, s->numFields);
3402: *subs = s->field[field];
3403: PetscFunctionReturn(PETSC_SUCCESS);
3404: }
3406: PetscClassId PETSC_SECTION_SYM_CLASSID;
3407: PetscFunctionList PetscSectionSymList = NULL;
3409: /*@
3410: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
3412: Collective
3414: Input Parameter:
3415: . comm - the MPI communicator
3417: Output Parameter:
3418: . sym - pointer to the new set of symmetries
3420: Level: developer
3422: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3423: @*/
3424: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3425: {
3426: PetscFunctionBegin;
3427: PetscAssertPointer(sym, 2);
3428: PetscCall(ISInitializePackage());
3430: PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3431: PetscFunctionReturn(PETSC_SUCCESS);
3432: }
3434: /*@
3435: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
3437: Collective
3439: Input Parameters:
3440: + sym - The section symmetry object
3441: - method - The name of the section symmetry type
3443: Level: developer
3445: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3446: @*/
3447: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3448: {
3449: PetscErrorCode (*r)(PetscSectionSym);
3450: PetscBool match;
3452: PetscFunctionBegin;
3454: PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3455: if (match) PetscFunctionReturn(PETSC_SUCCESS);
3457: PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3458: PetscCheck(r, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3459: PetscTryTypeMethod(sym, destroy);
3460: sym->ops->destroy = NULL;
3462: PetscCall((*r)(sym));
3463: PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3464: PetscFunctionReturn(PETSC_SUCCESS);
3465: }
3467: /*@
3468: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
3470: Not Collective
3472: Input Parameter:
3473: . sym - The section symmetry
3475: Output Parameter:
3476: . type - The index set type name
3478: Level: developer
3480: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3481: @*/
3482: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3483: {
3484: PetscFunctionBegin;
3486: PetscAssertPointer(type, 2);
3487: *type = ((PetscObject)sym)->type_name;
3488: PetscFunctionReturn(PETSC_SUCCESS);
3489: }
3491: /*@
3492: PetscSectionSymRegister - Registers a new section symmetry implementation
3494: Not Collective, No Fortran Support
3496: Input Parameters:
3497: + sname - The name of a new user-defined creation routine
3498: - function - The creation routine itself
3500: Level: developer
3502: Notes:
3503: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
3505: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3506: @*/
3507: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3508: {
3509: PetscFunctionBegin;
3510: PetscCall(ISInitializePackage());
3511: PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3512: PetscFunctionReturn(PETSC_SUCCESS);
3513: }
3515: /*@
3516: PetscSectionSymDestroy - Destroys a section symmetry.
3518: Collective
3520: Input Parameter:
3521: . sym - the section symmetry
3523: Level: developer
3525: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`
3526: @*/
3527: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3528: {
3529: SymWorkLink link, next;
3531: PetscFunctionBegin;
3532: if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3534: if (--((PetscObject)*sym)->refct > 0) {
3535: *sym = NULL;
3536: PetscFunctionReturn(PETSC_SUCCESS);
3537: }
3538: PetscTryTypeMethod(*sym, destroy);
3539: PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3540: for (link = (*sym)->workin; link; link = next) {
3541: PetscInt **perms = (PetscInt **)link->perms;
3542: PetscScalar **rots = (PetscScalar **)link->rots;
3543: PetscCall(PetscFree2(perms, rots));
3544: next = link->next;
3545: PetscCall(PetscFree(link));
3546: }
3547: (*sym)->workin = NULL;
3548: PetscCall(PetscHeaderDestroy(sym));
3549: PetscFunctionReturn(PETSC_SUCCESS);
3550: }
3552: /*@
3553: PetscSectionSymView - Displays a section symmetry
3555: Collective
3557: Input Parameters:
3558: + sym - the index set
3559: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3561: Level: developer
3563: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3564: @*/
3565: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3566: {
3567: PetscFunctionBegin;
3569: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3571: PetscCheckSameComm(sym, 1, viewer, 2);
3572: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3573: PetscTryTypeMethod(sym, view, viewer);
3574: PetscFunctionReturn(PETSC_SUCCESS);
3575: }
3577: /*@
3578: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3580: Collective
3582: Input Parameters:
3583: + section - the section describing data layout
3584: - sym - the symmetry describing the affect of orientation on the access of the data
3586: Level: developer
3588: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3589: @*/
3590: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3591: {
3592: PetscFunctionBegin;
3594: PetscCall(PetscSectionSymDestroy(§ion->sym));
3595: if (sym) {
3597: PetscCheckSameComm(section, 1, sym, 2);
3598: PetscCall(PetscObjectReference((PetscObject)sym));
3599: }
3600: section->sym = sym;
3601: PetscFunctionReturn(PETSC_SUCCESS);
3602: }
3604: /*@
3605: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3607: Not Collective
3609: Input Parameter:
3610: . section - the section describing data layout
3612: Output Parameter:
3613: . sym - the symmetry describing the affect of orientation on the access of the data, provided previously by `PetscSectionSetSym()`
3615: Level: developer
3617: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3618: @*/
3619: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3620: {
3621: PetscFunctionBegin;
3623: *sym = section->sym;
3624: PetscFunctionReturn(PETSC_SUCCESS);
3625: }
3627: /*@
3628: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3630: Collective
3632: Input Parameters:
3633: + section - the section describing data layout
3634: . field - the field number
3635: - sym - the symmetry describing the affect of orientation on the access of the data
3637: Level: developer
3639: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3640: @*/
3641: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3642: {
3643: PetscFunctionBegin;
3645: PetscSectionCheckValidField(field, section->numFields);
3646: PetscCall(PetscSectionSetSym(section->field[field], sym));
3647: PetscFunctionReturn(PETSC_SUCCESS);
3648: }
3650: /*@
3651: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3653: Collective
3655: Input Parameters:
3656: + section - the section describing data layout
3657: - field - the field number
3659: Output Parameter:
3660: . sym - the symmetry describing the affect of orientation on the access of the data
3662: Level: developer
3664: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3665: @*/
3666: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3667: {
3668: PetscFunctionBegin;
3670: PetscSectionCheckValidField(field, section->numFields);
3671: *sym = section->field[field]->sym;
3672: PetscFunctionReturn(PETSC_SUCCESS);
3673: }
3675: /*@
3676: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3678: Not Collective
3680: Input Parameters:
3681: + section - the section
3682: . numPoints - the number of points
3683: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3684: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3685: context, see `DMPlexGetConeOrientation()`).
3687: Output Parameters:
3688: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3689: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3690: identity).
3692: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3693: .vb
3694: const PetscInt **perms;
3695: const PetscScalar **rots;
3696: PetscInt lOffset;
3698: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3699: for (i = 0, lOffset = 0; i < numPoints; i++) {
3700: PetscInt point = points[2*i], dof, sOffset;
3701: const PetscInt *perm = perms ? perms[i] : NULL;
3702: const PetscScalar *rot = rots ? rots[i] : NULL;
3704: PetscSectionGetDof(section,point,&dof);
3705: PetscSectionGetOffset(section,point,&sOffset);
3707: if (perm) { for (j = 0; j < dof; j++) lArray[lOffset + perm[j]] = sArray[sOffset + j]; }
3708: else { for (j = 0; j < dof; j++) lArray[lOffset + j ] = sArray[sOffset + j]; }
3709: if (rot) { for (j = 0; j < dof; j++) lArray[lOffset + j ] *= rot[j]; }
3710: lOffset += dof;
3711: }
3712: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3713: .ve
3715: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3716: .vb
3717: const PetscInt **perms;
3718: const PetscScalar **rots;
3719: PetscInt lOffset;
3721: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3722: for (i = 0, lOffset = 0; i < numPoints; i++) {
3723: PetscInt point = points[2*i], dof, sOffset;
3724: const PetscInt *perm = perms ? perms[i] : NULL;
3725: const PetscScalar *rot = rots ? rots[i] : NULL;
3727: PetscSectionGetDof(section,point,&dof);
3728: PetscSectionGetOffset(section,point,&sOff);
3730: if (perm) { for (j = 0; j < dof; j++) sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.); }
3731: else { for (j = 0; j < dof; j++) sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.); }
3732: offset += dof;
3733: }
3734: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3735: .ve
3737: Level: developer
3739: Notes:
3740: `PetscSectionSetSym()` must have been previously called to provide the symmetries to the `PetscSection`
3742: Use `PetscSectionRestorePointSyms()` when finished with the data
3744: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3745: @*/
3746: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3747: {
3748: PetscSectionSym sym;
3750: PetscFunctionBegin;
3752: if (numPoints) PetscAssertPointer(points, 3);
3753: if (perms) *perms = NULL;
3754: if (rots) *rots = NULL;
3755: sym = section->sym;
3756: if (sym && (perms || rots)) {
3757: SymWorkLink link;
3759: if (sym->workin) {
3760: link = sym->workin;
3761: sym->workin = sym->workin->next;
3762: } else {
3763: PetscCall(PetscNew(&link));
3764: }
3765: if (numPoints > link->numPoints) {
3766: PetscInt **perms = (PetscInt **)link->perms;
3767: PetscScalar **rots = (PetscScalar **)link->rots;
3768: PetscCall(PetscFree2(perms, rots));
3769: PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3770: link->numPoints = numPoints;
3771: }
3772: link->next = sym->workout;
3773: sym->workout = link;
3774: PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3775: PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3776: PetscUseTypeMethod(sym, getpoints, section, numPoints, points, link->perms, link->rots);
3777: if (perms) *perms = link->perms;
3778: if (rots) *rots = link->rots;
3779: }
3780: PetscFunctionReturn(PETSC_SUCCESS);
3781: }
3783: /*@
3784: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3786: Not Collective
3788: Input Parameters:
3789: + section - the section
3790: . numPoints - the number of points
3791: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3792: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3793: context, see `DMPlexGetConeOrientation()`).
3794: . perms - The permutations for the given orientations: set to `NULL` at conclusion
3795: - rots - The field rotations symmetries for the given orientations: set to `NULL` at conclusion
3797: Level: developer
3799: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3800: @*/
3801: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3802: {
3803: PetscSectionSym sym;
3805: PetscFunctionBegin;
3807: sym = section->sym;
3808: if (sym && (perms || rots)) {
3809: SymWorkLink *p, link;
3811: for (p = &sym->workout; (link = *p); p = &link->next) {
3812: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3813: *p = link->next;
3814: link->next = sym->workin;
3815: sym->workin = link;
3816: if (perms) *perms = NULL;
3817: if (rots) *rots = NULL;
3818: PetscFunctionReturn(PETSC_SUCCESS);
3819: }
3820: }
3821: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3822: }
3823: PetscFunctionReturn(PETSC_SUCCESS);
3824: }
3826: /*@
3827: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3829: Not Collective
3831: Input Parameters:
3832: + section - the section
3833: . field - the field of the section
3834: . numPoints - the number of points
3835: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3836: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3837: context, see `DMPlexGetConeOrientation()`).
3839: Output Parameters:
3840: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3841: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3842: identity).
3844: Level: developer
3846: Notes:
3847: `PetscSectionSetFieldSym()` must have been previously called to provide the symmetries to the `PetscSection`
3849: Use `PetscSectionRestoreFieldPointSyms()` when finished with the data
3851: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3852: @*/
3853: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3854: {
3855: PetscFunctionBegin;
3857: PetscCheck(field <= section->numFields, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "field %" PetscInt_FMT " greater than number of fields (%" PetscInt_FMT ") in section", field, section->numFields);
3858: PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3859: PetscFunctionReturn(PETSC_SUCCESS);
3860: }
3862: /*@
3863: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3865: Not Collective
3867: Input Parameters:
3868: + section - the section
3869: . field - the field number
3870: . numPoints - the number of points
3871: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3872: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3873: context, see `DMPlexGetConeOrientation()`).
3874: . perms - The permutations for the given orientations: set to NULL at conclusion
3875: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3877: Level: developer
3879: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3880: @*/
3881: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3882: {
3883: PetscFunctionBegin;
3885: PetscCheck(field <= section->numFields, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "field %" PetscInt_FMT " greater than number of fields (%" PetscInt_FMT ") in section", field, section->numFields);
3886: PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3887: PetscFunctionReturn(PETSC_SUCCESS);
3888: }
3890: /*@
3891: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3893: Not Collective
3895: Input Parameter:
3896: . sym - the `PetscSectionSym`
3898: Output Parameter:
3899: . nsym - the equivalent symmetries
3901: Level: developer
3903: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3904: @*/
3905: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3906: {
3907: PetscFunctionBegin;
3910: PetscTryTypeMethod(sym, copy, nsym);
3911: PetscFunctionReturn(PETSC_SUCCESS);
3912: }
3914: /*@
3915: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3917: Collective
3919: Input Parameters:
3920: + sym - the `PetscSectionSym`
3921: - migrationSF - the distribution map from roots to leaves
3923: Output Parameter:
3924: . dsym - the redistributed symmetries
3926: Level: developer
3928: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3929: @*/
3930: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3931: {
3932: PetscFunctionBegin;
3935: PetscAssertPointer(dsym, 3);
3936: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3937: PetscFunctionReturn(PETSC_SUCCESS);
3938: }
3940: /*@
3941: PetscSectionGetUseFieldOffsets - Get the flag indicating if field offsets are used directly in a global section, rather than just the point offset
3943: Not Collective
3945: Input Parameter:
3946: . s - the global `PetscSection`
3948: Output Parameter:
3949: . flg - the flag
3951: Level: developer
3953: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3954: @*/
3955: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3956: {
3957: PetscFunctionBegin;
3959: *flg = s->useFieldOff;
3960: PetscFunctionReturn(PETSC_SUCCESS);
3961: }
3963: /*@
3964: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3966: Not Collective
3968: Input Parameters:
3969: + s - the global `PetscSection`
3970: - flg - the flag
3972: Level: developer
3974: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3975: @*/
3976: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3977: {
3978: PetscFunctionBegin;
3980: s->useFieldOff = flg;
3981: PetscFunctionReturn(PETSC_SUCCESS);
3982: }
3984: #define PetscSectionExpandPoints_Loop(TYPE) \
3985: do { \
3986: PetscInt i, n, o0, o1, size; \
3987: TYPE *a0 = (TYPE *)origArray, *a1; \
3988: PetscCall(PetscSectionGetStorageSize(s, &size)); \
3989: PetscCall(PetscMalloc1(size, &a1)); \
3990: for (i = 0; i < npoints; i++) { \
3991: PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
3992: PetscCall(PetscSectionGetOffset(s, i, &o1)); \
3993: PetscCall(PetscSectionGetDof(s, i, &n)); \
3994: PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n * unitsize)); \
3995: } \
3996: *newArray = (void *)a1; \
3997: } while (0)
3999: /*@
4000: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
4002: Not Collective
4004: Input Parameters:
4005: + origSection - the `PetscSection` describing the layout of the array
4006: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
4007: . origArray - the array; its size must be equal to the storage size of `origSection`
4008: - points - `IS` with points to extract; its indices must lie in the chart of `origSection`
4010: Output Parameters:
4011: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
4012: - newArray - the array of the extracted DOFs; its size is the storage size of `newSection`
4014: Level: developer
4016: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
4017: @*/
4018: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
4019: {
4020: PetscSection s;
4021: const PetscInt *points_;
4022: PetscInt i, n, npoints, pStart, pEnd;
4023: PetscMPIInt unitsize;
4025: PetscFunctionBegin;
4027: PetscAssertPointer(origArray, 3);
4029: if (newSection) PetscAssertPointer(newSection, 5);
4030: if (newArray) PetscAssertPointer(newArray, 6);
4031: PetscCallMPI(MPI_Type_size(dataType, &unitsize));
4032: PetscCall(ISGetLocalSize(points, &npoints));
4033: PetscCall(ISGetIndices(points, &points_));
4034: PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
4035: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
4036: PetscCall(PetscSectionSetChart(s, 0, npoints));
4037: for (i = 0; i < npoints; i++) {
4038: PetscCheck(points_[i] >= pStart && points_[i] < pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " (index %" PetscInt_FMT ") in input IS out of input section's chart", points_[i], i);
4039: PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
4040: PetscCall(PetscSectionSetDof(s, i, n));
4041: }
4042: PetscCall(PetscSectionSetUp(s));
4043: if (newArray) {
4044: if (dataType == MPIU_INT) {
4045: PetscSectionExpandPoints_Loop(PetscInt);
4046: } else if (dataType == MPIU_SCALAR) {
4047: PetscSectionExpandPoints_Loop(PetscScalar);
4048: } else if (dataType == MPIU_REAL) {
4049: PetscSectionExpandPoints_Loop(PetscReal);
4050: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
4051: }
4052: if (newSection) {
4053: *newSection = s;
4054: } else {
4055: PetscCall(PetscSectionDestroy(&s));
4056: }
4057: PetscCall(ISRestoreIndices(points, &points_));
4058: PetscFunctionReturn(PETSC_SUCCESS);
4059: }
4061: /*@
4062: PetscSectionMigrateData - Migrate data described by a `PetscSection` using a `PetscSF` that defines a original-to-new (root-to-leaf) point mapping
4064: Collective
4066: Input Parameters:
4067: + migratePointSF - defines the mapping (communication) of the root points to the leaf points
4068: . datatype - the type of data
4069: . rootSection - the `PetscSection` that describes the data layout on the root points (how many dof and what fields are associated with each root point)
4070: - rootData - the existing data array described by `rootSection`, may be `NULL` is storage size of `rootSection` is zero
4072: Output Parameters:
4073: + leafSection - the new `PetscSection` that describes the data layout on the leaf points
4074: . leafData - the redistributed data array that is associated with the leaf points
4075: - migrateDataSF - defines the mapping (communication) of the `rootData` array to the `leafData` array, may be `NULL` if not needed
4077: Level: advanced
4079: Notes:
4080: This function can best be thought of as applying `PetscSFBcastBegin()` to an array described by a `PetscSection`.
4081: While `PetscSFBcastBegin()` is limited to broadcasting data that is of the same size for every index, this function allows the data to be a different size for each index.
4082: The size and layout of that irregularly sized data before and after `PetscSFBcastBegin()` is described by the `rootSection` and `leafSection`, respectively.
4084: This function combines `PetscSFDistributeSection()`, `PetscSFCreateSectionSF()`, and `PetscSFBcastBegin()`/`PetscSFBcastEnd()` into a single call.
4085: `migrateDataSF` can be used to repeat the `PetscSFBcastBegin()`/`PetscSFBcastEnd()` on a different data array described by the same `rootSection`.
4087: This should not be used for global-to-local type communication patterns.
4088: For this use case, see `PetscSectionCreateGlobalSection()` and `PetscSFSetGraphSection()`.
4090: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSFDistributeSection()`, `PetscSFCreateSectionSF()`, `DMPlexDistributeData()`
4091: @*/
4092: PetscErrorCode PetscSectionMigrateData(PetscSF migratePointSF, MPI_Datatype datatype, PetscSection rootSection, const void *rootData, PetscSection leafSection, void *leafData[], PetscSF *migrateDataSF)
4093: {
4094: PetscSF fieldSF;
4095: PetscInt *remoteOffsets, fieldSize;
4096: PetscMPIInt dataSize;
4098: PetscFunctionBegin;
4101: if (rootData) PetscAssertPointer(rootData, 4);
4102: else {
4103: PetscInt size;
4104: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
4105: PetscCheck(size == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "originalData may be NULL iff the storage size of originalSection is zero, but is %" PetscInt_FMT, size);
4106: }
4108: PetscAssertPointer(leafData, 6);
4109: if (migrateDataSF) PetscAssertPointer(migrateDataSF, 7);
4111: PetscCall(PetscSFDistributeSection(migratePointSF, rootSection, &remoteOffsets, leafSection));
4112: PetscCall(PetscSFCreateSectionSF(migratePointSF, rootSection, remoteOffsets, leafSection, &fieldSF));
4113: PetscCall(PetscFree(remoteOffsets));
4115: PetscCall(PetscSectionGetStorageSize(leafSection, &fieldSize));
4116: PetscCallMPI(MPI_Type_size(datatype, &dataSize));
4117: PetscCall(PetscMalloc(fieldSize * dataSize, leafData));
4118: PetscCall(PetscSFBcastBegin(fieldSF, datatype, rootData, *leafData, MPI_REPLACE));
4119: PetscCall(PetscSFBcastEnd(fieldSF, datatype, rootData, *leafData, MPI_REPLACE));
4121: if (migrateDataSF) *migrateDataSF = fieldSF;
4122: else PetscCall(PetscSFDestroy(&fieldSF));
4123: PetscFunctionReturn(PETSC_SUCCESS);
4124: }