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](sec_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));
45: (*s)->pStart = -1;
46: (*s)->pEnd = -1;
47: (*s)->perm = NULL;
48: (*s)->pointMajor = PETSC_TRUE;
49: (*s)->includesConstraints = PETSC_TRUE;
50: (*s)->atlasDof = NULL;
51: (*s)->atlasOff = NULL;
52: (*s)->bc = NULL;
53: (*s)->bcIndices = NULL;
54: (*s)->setup = PETSC_FALSE;
55: (*s)->numFields = 0;
56: (*s)->fieldNames = NULL;
57: (*s)->field = NULL;
58: (*s)->useFieldOff = PETSC_FALSE;
59: (*s)->compNames = NULL;
60: (*s)->clObj = NULL;
61: (*s)->clHash = NULL;
62: (*s)->clSection = NULL;
63: (*s)->clPoints = NULL;
64: PetscCall(PetscSectionInvalidateMaxDof_Internal(*s));
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: /*@
69: PetscSectionCopy - Creates a shallow (if possible) copy of the `PetscSection`
71: Collective
73: Input Parameter:
74: . section - the `PetscSection`
76: Output Parameter:
77: . newSection - the copy
79: Level: intermediate
81: Developer Notes:
82: What exactly does shallow mean in this context?
84: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
85: @*/
86: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
87: {
88: PetscFunctionBegin;
91: PetscCall(PetscSectionCopy_Internal(section, newSection, NULL));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: PetscErrorCode PetscSectionCopy_Internal(PetscSection section, PetscSection newSection, PetscBT constrained_dofs)
96: {
97: PetscSectionSym sym;
98: IS perm;
99: PetscInt numFields, f, c, pStart, pEnd, p;
101: PetscFunctionBegin;
104: PetscCall(PetscSectionReset(newSection));
105: PetscCall(PetscSectionGetNumFields(section, &numFields));
106: if (numFields) PetscCall(PetscSectionSetNumFields(newSection, numFields));
107: for (f = 0; f < numFields; ++f) {
108: const char *fieldName = NULL, *compName = NULL;
109: PetscInt numComp = 0;
111: PetscCall(PetscSectionGetFieldName(section, f, &fieldName));
112: PetscCall(PetscSectionSetFieldName(newSection, f, fieldName));
113: PetscCall(PetscSectionGetFieldComponents(section, f, &numComp));
114: PetscCall(PetscSectionSetFieldComponents(newSection, f, numComp));
115: for (c = 0; c < numComp; ++c) {
116: PetscCall(PetscSectionGetComponentName(section, f, c, &compName));
117: PetscCall(PetscSectionSetComponentName(newSection, f, c, compName));
118: }
119: PetscCall(PetscSectionGetFieldSym(section, f, &sym));
120: PetscCall(PetscSectionSetFieldSym(newSection, f, sym));
121: }
122: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
123: PetscCall(PetscSectionSetChart(newSection, pStart, pEnd));
124: PetscCall(PetscSectionGetPermutation(section, &perm));
125: PetscCall(PetscSectionSetPermutation(newSection, perm));
126: PetscCall(PetscSectionGetSym(section, &sym));
127: PetscCall(PetscSectionSetSym(newSection, sym));
128: for (p = pStart; p < pEnd; ++p) {
129: PetscInt dof, cdof, fcdof = 0;
130: PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));
132: PetscCall(PetscSectionGetDof(section, p, &dof));
133: PetscCall(PetscSectionSetDof(newSection, p, dof));
134: if (force_constrained) cdof = dof;
135: else PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
136: if (cdof) PetscCall(PetscSectionSetConstraintDof(newSection, p, cdof));
137: for (f = 0; f < numFields; ++f) {
138: PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
139: PetscCall(PetscSectionSetFieldDof(newSection, p, f, dof));
140: if (cdof) {
141: if (force_constrained) fcdof = dof;
142: else PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
143: if (fcdof) PetscCall(PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof));
144: }
145: }
146: }
147: PetscCall(PetscSectionSetUp(newSection));
148: for (p = pStart; p < pEnd; ++p) {
149: PetscInt off, cdof, fcdof = 0;
150: const PetscInt *cInd;
151: PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));
153: /* Must set offsets in case they do not agree with the prefix sums */
154: PetscCall(PetscSectionGetOffset(section, p, &off));
155: PetscCall(PetscSectionSetOffset(newSection, p, off));
156: PetscCall(PetscSectionGetConstraintDof(newSection, p, &cdof));
157: if (cdof) {
158: if (force_constrained) cInd = NULL;
159: else PetscCall(PetscSectionGetConstraintIndices(section, p, &cInd));
160: PetscCall(PetscSectionSetConstraintIndices(newSection, p, cInd));
161: for (f = 0; f < numFields; ++f) {
162: PetscCall(PetscSectionGetFieldOffset(section, p, f, &off));
163: PetscCall(PetscSectionSetFieldOffset(newSection, p, f, off));
164: PetscCall(PetscSectionGetFieldConstraintDof(newSection, p, f, &fcdof));
165: if (fcdof) {
166: if (force_constrained) cInd = NULL;
167: else PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &cInd));
168: PetscCall(PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd));
169: }
170: }
171: }
172: }
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*@
177: PetscSectionClone - Creates a shallow (if possible) copy of the `PetscSection`
179: Collective
181: Input Parameter:
182: . section - the `PetscSection`
184: Output Parameter:
185: . newSection - the copy
187: Level: beginner
189: Developer Notes:
190: With standard PETSc terminology this should be called `PetscSectionDuplicate()`
192: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
193: @*/
194: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
195: {
196: PetscFunctionBegin;
198: PetscAssertPointer(newSection, 2);
199: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection));
200: PetscCall(PetscSectionCopy(section, *newSection));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*@
205: PetscSectionSetFromOptions - sets parameters in a `PetscSection` from the options database
207: Collective
209: Input Parameter:
210: . s - the `PetscSection`
212: Options Database Key:
213: . -petscsection_point_major - `PETSC_TRUE` for point-major order
215: Level: intermediate
217: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
218: @*/
219: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
220: {
221: PetscFunctionBegin;
223: PetscObjectOptionsBegin((PetscObject)s);
224: PetscCall(PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL));
225: /* process any options handlers added with PetscObjectAddOptionsHandler() */
226: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject));
227: PetscOptionsEnd();
228: PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view"));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*@
233: PetscSectionCompare - Compares two sections
235: Collective
237: Input Parameters:
238: + s1 - the first `PetscSection`
239: - s2 - the second `PetscSection`
241: Output Parameter:
242: . congruent - `PETSC_TRUE` if the two sections are congruent, `PETSC_FALSE` otherwise
244: Level: intermediate
246: Note:
247: Field names are disregarded.
249: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
250: @*/
251: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
252: {
253: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
254: const PetscInt *idx1, *idx2;
255: IS perm1, perm2;
256: PetscBool flg;
257: PetscMPIInt mflg;
259: PetscFunctionBegin;
262: PetscAssertPointer(congruent, 3);
263: flg = PETSC_FALSE;
265: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg));
266: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
267: *congruent = PETSC_FALSE;
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: PetscCall(PetscSectionGetChart(s1, &pStart, &pEnd));
272: PetscCall(PetscSectionGetChart(s2, &n1, &n2));
273: if (pStart != n1 || pEnd != n2) goto not_congruent;
275: PetscCall(PetscSectionGetPermutation(s1, &perm1));
276: PetscCall(PetscSectionGetPermutation(s2, &perm2));
277: if (perm1 && perm2) {
278: PetscCall(ISEqual(perm1, perm2, congruent));
279: if (!*congruent) goto not_congruent;
280: } else if (perm1 != perm2) goto not_congruent;
282: for (p = pStart; p < pEnd; ++p) {
283: PetscCall(PetscSectionGetOffset(s1, p, &n1));
284: PetscCall(PetscSectionGetOffset(s2, p, &n2));
285: if (n1 != n2) goto not_congruent;
287: PetscCall(PetscSectionGetDof(s1, p, &n1));
288: PetscCall(PetscSectionGetDof(s2, p, &n2));
289: if (n1 != n2) goto not_congruent;
291: PetscCall(PetscSectionGetConstraintDof(s1, p, &ncdof));
292: PetscCall(PetscSectionGetConstraintDof(s2, p, &n2));
293: if (ncdof != n2) goto not_congruent;
295: PetscCall(PetscSectionGetConstraintIndices(s1, p, &idx1));
296: PetscCall(PetscSectionGetConstraintIndices(s2, p, &idx2));
297: PetscCall(PetscArraycmp(idx1, idx2, ncdof, congruent));
298: if (!*congruent) goto not_congruent;
299: }
301: PetscCall(PetscSectionGetNumFields(s1, &nfields));
302: PetscCall(PetscSectionGetNumFields(s2, &n2));
303: if (nfields != n2) goto not_congruent;
305: for (f = 0; f < nfields; ++f) {
306: PetscCall(PetscSectionGetFieldComponents(s1, f, &n1));
307: PetscCall(PetscSectionGetFieldComponents(s2, f, &n2));
308: if (n1 != n2) goto not_congruent;
310: for (p = pStart; p < pEnd; ++p) {
311: PetscCall(PetscSectionGetFieldOffset(s1, p, f, &n1));
312: PetscCall(PetscSectionGetFieldOffset(s2, p, f, &n2));
313: if (n1 != n2) goto not_congruent;
315: PetscCall(PetscSectionGetFieldDof(s1, p, f, &n1));
316: PetscCall(PetscSectionGetFieldDof(s2, p, f, &n2));
317: if (n1 != n2) goto not_congruent;
319: PetscCall(PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof));
320: PetscCall(PetscSectionGetFieldConstraintDof(s2, p, f, &n2));
321: if (nfcdof != n2) goto not_congruent;
323: PetscCall(PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1));
324: PetscCall(PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2));
325: PetscCall(PetscArraycmp(idx1, idx2, nfcdof, congruent));
326: if (!*congruent) goto not_congruent;
327: }
328: }
330: flg = PETSC_TRUE;
331: not_congruent:
332: PetscCall(MPIU_Allreduce(&flg, congruent, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1)));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /*@
337: PetscSectionGetNumFields - Returns the number of fields in a `PetscSection`, or 0 if no fields were defined.
339: Not Collective
341: Input Parameter:
342: . s - the `PetscSection`
344: Output Parameter:
345: . numFields - the number of fields defined, or 0 if none were defined
347: Level: intermediate
349: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
350: @*/
351: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
352: {
353: PetscFunctionBegin;
355: PetscAssertPointer(numFields, 2);
356: *numFields = s->numFields;
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: /*@
361: PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`
363: Not Collective
365: Input Parameters:
366: + s - the `PetscSection`
367: - numFields - the number of fields
369: Level: intermediate
371: Notes:
372: Calling this destroys all the information in the `PetscSection` including the chart.
374: You must call `PetscSectionSetChart()` after calling this.
376: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetNumFields()`, `PetscSectionSetChart()`, `PetscSectionReset()`
377: @*/
378: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
379: {
380: PetscInt f;
382: PetscFunctionBegin;
384: PetscCheck(numFields > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %" PetscInt_FMT " must be positive", numFields);
385: PetscCall(PetscSectionReset(s));
387: s->numFields = numFields;
388: PetscCall(PetscMalloc1(s->numFields, &s->numFieldComponents));
389: PetscCall(PetscMalloc1(s->numFields, &s->fieldNames));
390: PetscCall(PetscMalloc1(s->numFields, &s->compNames));
391: PetscCall(PetscMalloc1(s->numFields, &s->field));
392: for (f = 0; f < s->numFields; ++f) {
393: char name[64];
395: s->numFieldComponents[f] = 1;
397: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]));
398: PetscCall(PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f));
399: PetscCall(PetscStrallocpy(name, (char **)&s->fieldNames[f]));
400: PetscCall(PetscSNPrintf(name, 64, "Component_0"));
401: PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]));
402: PetscCall(PetscStrallocpy(name, (char **)&s->compNames[f][0]));
403: }
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: /*@C
408: PetscSectionGetFieldName - Returns the name of a field in the `PetscSection`
410: Not Collective
412: Input Parameters:
413: + s - the `PetscSection`
414: - field - the field number
416: Output Parameter:
417: . fieldName - the field name
419: Level: intermediate
421: Note:
422: Will error if the field number is out of range
424: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
425: @*/
426: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
427: {
428: PetscFunctionBegin;
430: PetscAssertPointer(fieldName, 3);
431: PetscSectionCheckValidField(field, s->numFields);
432: *fieldName = s->fieldNames[field];
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: /*@C
437: PetscSectionSetFieldName - Sets the name of a field in the `PetscSection`
439: Not Collective
441: Input Parameters:
442: + s - the `PetscSection`
443: . field - the field number
444: - fieldName - the field name
446: Level: intermediate
448: Note:
449: Will error if the field number is out of range
451: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
452: @*/
453: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
454: {
455: PetscFunctionBegin;
457: if (fieldName) PetscAssertPointer(fieldName, 3);
458: PetscSectionCheckValidField(field, s->numFields);
459: PetscCall(PetscFree(s->fieldNames[field]));
460: PetscCall(PetscStrallocpy(fieldName, (char **)&s->fieldNames[field]));
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@C
465: PetscSectionGetComponentName - Gets the name of a field component in the `PetscSection`
467: Not Collective
469: Input Parameters:
470: + s - the `PetscSection`
471: . field - the field number
472: - comp - the component number
474: Output Parameter:
475: . compName - the component name
477: Level: intermediate
479: Note:
480: Will error if the field or component number do not exist
482: Developer Notes:
483: The function name should have Field in it since they are field components.
485: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
486: `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
487: @*/
488: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
489: {
490: PetscFunctionBegin;
492: PetscAssertPointer(compName, 4);
493: PetscSectionCheckValidField(field, s->numFields);
494: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
495: *compName = s->compNames[field][comp];
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /*@C
500: PetscSectionSetComponentName - Sets the name of a field component in the `PetscSection`
502: Not Collective
504: Input Parameters:
505: + s - the `PetscSection`
506: . field - the field number
507: . comp - the component number
508: - compName - the component name
510: Level: advanced
512: Note:
513: Will error if the field or component number do not exist
515: Developer Notes:
516: The function name should have Field in it since they are field components.
518: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
519: `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
520: @*/
521: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
522: {
523: PetscFunctionBegin;
525: if (compName) PetscAssertPointer(compName, 4);
526: PetscSectionCheckValidField(field, s->numFields);
527: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
528: PetscCall(PetscFree(s->compNames[field][comp]));
529: PetscCall(PetscStrallocpy(compName, (char **)&s->compNames[field][comp]));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /*@
534: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
536: Not Collective
538: Input Parameters:
539: + s - the `PetscSection`
540: - field - the field number
542: Output Parameter:
543: . numComp - the number of field components
545: Level: advanced
547: Developer Notes:
548: This function is misnamed. There is a Num in `PetscSectionGetNumFields()` but not in this name
550: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`,
551: `PetscSectionSetComponentName()`, `PetscSectionGetComponentName()`
552: @*/
553: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
554: {
555: PetscFunctionBegin;
557: PetscAssertPointer(numComp, 3);
558: PetscSectionCheckValidField(field, s->numFields);
559: *numComp = s->numFieldComponents[field];
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: /*@
564: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
566: Not Collective
568: Input Parameters:
569: + s - the `PetscSection`
570: . field - the field number
571: - numComp - the number of field components
573: Level: advanced
575: Note:
576: This number can be different than the values set with `PetscSectionSetFieldDof()`. It can be used to indicate the number of
577: components in the field of the underlying physical model which may be different than the number of degrees of freedom needed
578: 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
579: an face based model for velocity (where the velocity normal to the face is stored) there is only 1 dof for each face point.
581: The value set with this function are not needed or used in `PetscSectionSetUp()`.
583: Developer Notes:
584: This function is misnamed. There is a Num in `PetscSectionSetNumFields()` but not in this name
586: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionSetComponentName()`,
587: `PetscSectionGetComponentName()`, `PetscSectionGetNumFields()`
588: @*/
589: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
590: {
591: PetscInt c;
593: PetscFunctionBegin;
595: PetscSectionCheckValidField(field, s->numFields);
596: if (s->compNames) {
597: for (c = 0; c < s->numFieldComponents[field]; ++c) PetscCall(PetscFree(s->compNames[field][c]));
598: PetscCall(PetscFree(s->compNames[field]));
599: }
601: s->numFieldComponents[field] = numComp;
602: if (numComp) {
603: PetscCall(PetscMalloc1(numComp, (char ***)&s->compNames[field]));
604: for (c = 0; c < numComp; ++c) {
605: char name[64];
607: PetscCall(PetscSNPrintf(name, 64, "%" PetscInt_FMT, c));
608: PetscCall(PetscStrallocpy(name, (char **)&s->compNames[field][c]));
609: }
610: }
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: /*@
615: PetscSectionGetChart - Returns the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
617: Not Collective
619: Input Parameter:
620: . s - the `PetscSection`
622: Output Parameters:
623: + pStart - the first point
624: - pEnd - one past the last point
626: Level: intermediate
628: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
629: @*/
630: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
631: {
632: PetscFunctionBegin;
634: if (pStart) *pStart = s->pStart;
635: if (pEnd) *pEnd = s->pEnd;
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: /*@
640: PetscSectionSetChart - Sets the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
642: Not Collective
644: Input Parameters:
645: + s - the `PetscSection`
646: . pStart - the first point
647: - pEnd - one past the last point, `pStart` $ \le $ `pEnd`
649: Level: intermediate
651: Notes:
652: The charts on different MPI processes may (and often do) overlap
654: If you intend to use `PetscSectionSetNumFields()` it must be called before this call.
656: The chart for all fields created with `PetscSectionSetNumFields()` is the same as this chart.
658: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`, `PetscSectionSetNumFields()`
659: @*/
660: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
661: {
662: PetscInt f;
664: PetscFunctionBegin;
666: PetscCheck(pEnd >= pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Chart pEnd %" PetscInt_FMT " cannot be smaller than chart pStart %" PetscInt_FMT, pEnd, pStart);
667: if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(PETSC_SUCCESS);
668: /* Cannot Reset() because it destroys field information */
669: s->setup = PETSC_FALSE;
670: PetscCall(PetscSectionDestroy(&s->bc));
671: PetscCall(PetscFree(s->bcIndices));
672: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
674: s->pStart = pStart;
675: s->pEnd = pEnd;
676: PetscCall(PetscMalloc2((pEnd - pStart), &s->atlasDof, (pEnd - pStart), &s->atlasOff));
677: PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart));
678: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: /*@
683: PetscSectionGetPermutation - Returns the permutation of [0, `pEnd` - `pStart`) or `NULL` that was set with `PetscSectionSetPermutation()`
685: Not Collective
687: Input Parameter:
688: . s - the `PetscSection`
690: Output Parameter:
691: . perm - The permutation as an `IS`
693: Level: intermediate
695: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
696: @*/
697: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
698: {
699: PetscFunctionBegin;
701: if (perm) {
702: PetscAssertPointer(perm, 2);
703: *perm = s->perm;
704: }
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: /*@
709: PetscSectionSetPermutation - Sets a permutation of the chart for this section, [0, `pEnd` - `pStart`), which determines the order to store the `PetscSection` information
711: Not Collective
713: Input Parameters:
714: + s - the `PetscSection`
715: - perm - the permutation of points
717: Level: intermediate
719: Notes:
720: The permutation must be provided before `PetscSectionSetUp()`.
722: The data in the `PetscSection` are permuted but the access via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` is not changed
724: Compare to `PetscSectionPermute()`
726: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetUp()`, `PetscSectionGetPermutation()`, `PetscSectionPermute()`, `PetscSectionCreate()`
727: @*/
728: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
729: {
730: PetscFunctionBegin;
733: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
734: if (s->perm != perm) {
735: PetscCall(ISDestroy(&s->perm));
736: if (perm) {
737: s->perm = perm;
738: PetscCall(PetscObjectReference((PetscObject)s->perm));
739: }
740: }
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: /*@C
745: PetscSectionGetBlockStarts - Returns a table indicating which points start new blocks
747: Not Collective
749: Input Parameter:
750: . s - the `PetscSection`
752: Output Parameter:
753: . blockStarts - The `PetscBT` with a 1 for each point that begins a block
755: Notes:
756: The table is on [0, `pEnd` - `pStart`).
758: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
760: Level: intermediate
762: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
763: @*/
764: PetscErrorCode PetscSectionGetBlockStarts(PetscSection s, PetscBT *blockStarts)
765: {
766: PetscFunctionBegin;
768: if (blockStarts) {
769: PetscAssertPointer(blockStarts, 2);
770: *blockStarts = s->blockStarts;
771: }
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: /*@C
776: PetscSectionSetBlockStarts - Sets a table indicating which points start new blocks
778: Not Collective
780: Input Parameters:
781: + s - the `PetscSection`
782: - blockStarts - The `PetscBT` with a 1 for each point that begins a block
784: Level: intermediate
786: Notes:
787: 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.
789: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
791: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionGetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
792: @*/
793: PetscErrorCode PetscSectionSetBlockStarts(PetscSection s, PetscBT blockStarts)
794: {
795: PetscFunctionBegin;
797: if (s->blockStarts != blockStarts) {
798: PetscCall(PetscBTDestroy(&s->blockStarts));
799: s->blockStarts = blockStarts;
800: }
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: /*@
805: PetscSectionGetPointMajor - Returns the flag for dof ordering, `PETSC_TRUE` if it is point major, `PETSC_FALSE` if it is field major
807: Not Collective
809: Input Parameter:
810: . s - the `PetscSection`
812: Output Parameter:
813: . pm - the flag for point major ordering
815: Level: intermediate
817: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetPointMajor()`
818: @*/
819: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
820: {
821: PetscFunctionBegin;
823: PetscAssertPointer(pm, 2);
824: *pm = s->pointMajor;
825: PetscFunctionReturn(PETSC_SUCCESS);
826: }
828: /*@
829: PetscSectionSetPointMajor - Sets the flag for dof ordering, `PETSC_TRUE` for point major, otherwise it will be field major
831: Not Collective
833: Input Parameters:
834: + s - the `PetscSection`
835: - pm - the flag for point major ordering
837: Level: intermediate
839: Note:
840: 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.
842: Point major order means the degrees of freedom are stored as follows
843: .vb
844: all the degrees of freedom for each point are stored contiguously, one point after another (respecting a permutation set with `PetscSectionSetPermutation()`)
845: for each point
846: the degrees of freedom for each field (starting with the unnamed default field) are listed in order by field
847: .ve
849: Field major order means the degrees of freedom are stored as follows
850: .vb
851: all degrees of freedom for each field (including the unnamed default field) are stored contiguously, one field after another
852: for each field (started with unnamed default field)
853: the degrees of freedom for each point are listed in order by point (respecting a permutation set with `PetscSectionSetPermutation()`)
854: .ve
856: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`, `PetscSectionSetPermutation()`
857: @*/
858: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
859: {
860: PetscFunctionBegin;
862: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
863: s->pointMajor = pm;
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets in the `PetscSection`.
869: The value is set with `PetscSectionSetIncludesConstraints()`
871: Not Collective
873: Input Parameter:
874: . s - the `PetscSection`
876: Output Parameter:
877: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
879: Level: intermediate
881: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
882: @*/
883: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
884: {
885: PetscFunctionBegin;
887: PetscAssertPointer(includesConstraints, 2);
888: *includesConstraints = s->includesConstraints;
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: /*@
893: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
895: Not Collective
897: Input Parameters:
898: + s - the `PetscSection`
899: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
901: Level: intermediate
903: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
904: @*/
905: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
906: {
907: PetscFunctionBegin;
909: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
910: s->includesConstraints = includesConstraints;
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: /*@
915: PetscSectionGetDof - Return the total number of degrees of freedom associated with a given point.
917: Not Collective
919: Input Parameters:
920: + s - the `PetscSection`
921: - point - the point
923: Output Parameter:
924: . numDof - the number of dof
926: Level: intermediate
928: Notes:
929: In a global section, this size will be negative for points not owned by this process.
931: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
933: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
934: @*/
935: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
936: {
937: PetscFunctionBeginHot;
939: PetscAssertPointer(numDof, 3);
940: 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);
941: *numDof = s->atlasDof[point - s->pStart];
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: /*@
946: PetscSectionSetDof - Sets the total number of degrees of freedom associated with a given point.
948: Not Collective
950: Input Parameters:
951: + s - the `PetscSection`
952: . point - the point
953: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
955: Level: intermediate
957: Note:
958: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
960: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
961: @*/
962: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
963: {
964: PetscFunctionBegin;
966: 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);
967: s->atlasDof[point - s->pStart] = numDof;
968: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
969: PetscFunctionReturn(PETSC_SUCCESS);
970: }
972: /*@
973: PetscSectionAddDof - Adds to the total number of degrees of freedom associated with a given point.
975: Not Collective
977: Input Parameters:
978: + s - the `PetscSection`
979: . point - the point
980: - numDof - the number of additional dof
982: Level: intermediate
984: Note:
985: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
987: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
988: @*/
989: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
990: {
991: PetscFunctionBeginHot;
993: 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);
994: PetscCheck(numDof >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numDof %" PetscInt_FMT " should not be negative", numDof);
995: s->atlasDof[point - s->pStart] += numDof;
996: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /*@
1001: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
1003: Not Collective
1005: Input Parameters:
1006: + s - the `PetscSection`
1007: . point - the point
1008: - field - the field
1010: Output Parameter:
1011: . numDof - the number of dof
1013: Level: intermediate
1015: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
1016: @*/
1017: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1018: {
1019: PetscFunctionBegin;
1021: PetscAssertPointer(numDof, 4);
1022: PetscSectionCheckValidField(field, s->numFields);
1023: PetscCall(PetscSectionGetDof(s->field[field], point, numDof));
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }
1027: /*@
1028: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
1030: Not Collective
1032: Input Parameters:
1033: + s - the `PetscSection`
1034: . point - the point
1035: . field - the field
1036: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
1038: Level: intermediate
1040: Note:
1041: 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
1042: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1044: This is equivalent to
1045: .vb
1046: PetscSection fs;
1047: PetscSectionGetField(s,field,&fs)
1048: PetscSectionSetDof(fs,numDof)
1049: .ve
1051: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`, `PetscSectionAddDof()`, `PetscSectionSetDof()`
1052: @*/
1053: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1054: {
1055: PetscFunctionBegin;
1057: PetscSectionCheckValidField(field, s->numFields);
1058: PetscCall(PetscSectionSetDof(s->field[field], point, numDof));
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: /*@
1063: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
1065: Not Collective
1067: Input Parameters:
1068: + s - the `PetscSection`
1069: . point - the point
1070: . field - the field
1071: - numDof - the number of dof
1073: Level: intermediate
1075: Notes:
1076: 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
1077: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1079: This is equivalent to
1080: .vb
1081: PetscSection fs;
1082: PetscSectionGetField(s,field,&fs)
1083: PetscSectionAddDof(fs,numDof)
1084: .ve
1086: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
1087: @*/
1088: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1089: {
1090: PetscFunctionBegin;
1092: PetscSectionCheckValidField(field, s->numFields);
1093: PetscCall(PetscSectionAddDof(s->field[field], point, numDof));
1094: PetscFunctionReturn(PETSC_SUCCESS);
1095: }
1097: /*@
1098: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
1100: Not Collective
1102: Input Parameters:
1103: + s - the `PetscSection`
1104: - point - the point
1106: Output Parameter:
1107: . numDof - the number of dof which are fixed by constraints
1109: Level: intermediate
1111: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
1112: @*/
1113: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
1114: {
1115: PetscFunctionBegin;
1117: PetscAssertPointer(numDof, 3);
1118: if (s->bc) {
1119: PetscCall(PetscSectionGetDof(s->bc, point, numDof));
1120: } else *numDof = 0;
1121: PetscFunctionReturn(PETSC_SUCCESS);
1122: }
1124: /*@
1125: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
1127: Not Collective
1129: Input Parameters:
1130: + s - the `PetscSection`
1131: . point - the point
1132: - numDof - the number of dof which are fixed by constraints
1134: Level: intermediate
1136: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1137: @*/
1138: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1139: {
1140: PetscFunctionBegin;
1142: if (numDof) {
1143: PetscCall(PetscSectionCheckConstraints_Private(s));
1144: PetscCall(PetscSectionSetDof(s->bc, point, numDof));
1145: }
1146: PetscFunctionReturn(PETSC_SUCCESS);
1147: }
1149: /*@
1150: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
1152: Not Collective
1154: Input Parameters:
1155: + s - the `PetscSection`
1156: . point - the point
1157: - numDof - the number of additional dof which are fixed by constraints
1159: Level: intermediate
1161: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1162: @*/
1163: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1164: {
1165: PetscFunctionBegin;
1167: if (numDof) {
1168: PetscCall(PetscSectionCheckConstraints_Private(s));
1169: PetscCall(PetscSectionAddDof(s->bc, point, numDof));
1170: }
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: /*@
1175: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
1177: Not Collective
1179: Input Parameters:
1180: + s - the `PetscSection`
1181: . point - the point
1182: - field - the field
1184: Output Parameter:
1185: . numDof - the number of dof which are fixed by constraints
1187: Level: intermediate
1189: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1190: @*/
1191: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1192: {
1193: PetscFunctionBegin;
1195: PetscAssertPointer(numDof, 4);
1196: PetscSectionCheckValidField(field, s->numFields);
1197: PetscCall(PetscSectionGetConstraintDof(s->field[field], point, numDof));
1198: PetscFunctionReturn(PETSC_SUCCESS);
1199: }
1201: /*@
1202: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1204: Not Collective
1206: Input Parameters:
1207: + s - the `PetscSection`
1208: . point - the point
1209: . field - the field
1210: - numDof - the number of dof which are fixed by constraints
1212: Level: intermediate
1214: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1215: @*/
1216: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1217: {
1218: PetscFunctionBegin;
1220: PetscSectionCheckValidField(field, s->numFields);
1221: PetscCall(PetscSectionSetConstraintDof(s->field[field], point, numDof));
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: /*@
1226: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1228: Not Collective
1230: Input Parameters:
1231: + s - the `PetscSection`
1232: . point - the point
1233: . field - the field
1234: - numDof - the number of additional dof which are fixed by constraints
1236: Level: intermediate
1238: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1239: @*/
1240: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1241: {
1242: PetscFunctionBegin;
1244: PetscSectionCheckValidField(field, s->numFields);
1245: PetscCall(PetscSectionAddConstraintDof(s->field[field], point, numDof));
1246: PetscFunctionReturn(PETSC_SUCCESS);
1247: }
1249: /*@
1250: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1252: Not Collective
1254: Input Parameter:
1255: . s - the `PetscSection`
1257: Level: advanced
1259: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1260: @*/
1261: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1262: {
1263: PetscFunctionBegin;
1265: if (s->bc) {
1266: const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;
1268: PetscCall(PetscSectionSetUp(s->bc));
1269: PetscCall(PetscMalloc1((last >= 0 ? s->bc->atlasOff[last] + s->bc->atlasDof[last] : 0), &s->bcIndices));
1270: }
1271: PetscFunctionReturn(PETSC_SUCCESS);
1272: }
1274: /*@
1275: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point in preparation for use of the `PetscSection`
1277: Not Collective
1279: Input Parameter:
1280: . s - the `PetscSection`
1282: Level: intermediate
1284: Notes:
1285: If used, `PetscSectionSetPermutation()` must be called before this routine.
1287: `PetscSectionSetPointMajor()`, cannot be called after this routine.
1289: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
1290: @*/
1291: PetscErrorCode PetscSectionSetUp(PetscSection s)
1292: {
1293: const PetscInt *pind = NULL;
1294: PetscInt offset = 0, foff, p, f;
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: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1306: const PetscInt q = pind ? pind[p] : p;
1308: /* Set point offset */
1309: s->atlasOff[q] = offset;
1310: offset += s->atlasDof[q];
1311: /* Set field offset */
1312: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1313: PetscSection sf = s->field[f];
1315: sf->atlasOff[q] = foff;
1316: foff += sf->atlasDof[q];
1317: }
1318: }
1319: } else {
1320: /* Set field offsets for all points */
1321: for (f = 0; f < s->numFields; ++f) {
1322: PetscSection sf = s->field[f];
1324: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1325: const PetscInt q = pind ? pind[p] : p;
1327: sf->atlasOff[q] = offset;
1328: offset += sf->atlasDof[q];
1329: }
1330: }
1331: /* Disable point offsets since these are unused */
1332: for (p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1333: }
1334: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1335: /* Setup BC sections */
1336: PetscCall(PetscSectionSetUpBC(s));
1337: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetUpBC(s->field[f]));
1338: PetscFunctionReturn(PETSC_SUCCESS);
1339: }
1341: /*@
1342: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the `PetscSection`
1344: Not Collective
1346: Input Parameter:
1347: . s - the `PetscSection`
1349: Output Parameter:
1350: . maxDof - the maximum dof
1352: Level: intermediate
1354: Notes:
1355: The returned number is up-to-date without need for `PetscSectionSetUp()`.
1357: 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
1358: the maximum over all points of the value returned by `PetscSectionGetDof()` on this MPI process
1360: Developer Notes:
1361: The returned number is calculated lazily and stashed.
1363: A call to `PetscSectionInvalidateMaxDof_Internal()` invalidates the stashed value.
1365: `PetscSectionInvalidateMaxDof_Internal()` is called in `PetscSectionSetDof()`, `PetscSectionAddDof()` and `PetscSectionReset()`
1367: It should also be called every time `atlasDof` is modified directly.
1369: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1370: @*/
1371: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1372: {
1373: PetscInt p;
1375: PetscFunctionBegin;
1377: PetscAssertPointer(maxDof, 2);
1378: if (s->maxDof == PETSC_MIN_INT) {
1379: s->maxDof = 0;
1380: for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1381: }
1382: *maxDof = s->maxDof;
1383: PetscFunctionReturn(PETSC_SUCCESS);
1384: }
1386: /*@
1387: PetscSectionGetStorageSize - Return the size of an array or local `Vec` capable of holding all the degrees of freedom defined in a `PetscSection`
1389: Not Collective
1391: Input Parameter:
1392: . s - the `PetscSection`
1394: Output Parameter:
1395: . size - the size of an array which can hold all the dofs
1397: Level: intermediate
1399: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1400: @*/
1401: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1402: {
1403: PetscInt p, n = 0;
1405: PetscFunctionBegin;
1407: PetscAssertPointer(size, 2);
1408: for (p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1409: *size = n;
1410: PetscFunctionReturn(PETSC_SUCCESS);
1411: }
1413: /*@
1414: PetscSectionGetConstrainedStorageSize - Return the size of an array or local `Vec` capable of holding all unconstrained degrees of freedom in a `PetscSection`
1416: Not Collective
1418: Input Parameter:
1419: . s - the `PetscSection`
1421: Output Parameter:
1422: . size - the size of an array which can hold all unconstrained dofs
1424: Level: intermediate
1426: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1427: @*/
1428: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1429: {
1430: PetscInt p, n = 0;
1432: PetscFunctionBegin;
1434: PetscAssertPointer(size, 2);
1435: for (p = 0; p < s->pEnd - s->pStart; ++p) {
1436: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1437: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1438: }
1439: *size = n;
1440: PetscFunctionReturn(PETSC_SUCCESS);
1441: }
1443: /*@
1444: PetscSectionCreateGlobalSection - Create a parallel section describing the global layout using
1445: a local (sequential) `PetscSection` on each MPI process and a `PetscSF` describing the section point overlap.
1447: Input Parameters:
1448: + s - The `PetscSection` for the local field layout
1449: . sf - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1450: . usePermutation - By default this is `PETSC_TRUE`, meaning any permutation of the local section is transferred to the global section
1451: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1452: - localOffsets - If `PETSC_TRUE`, use local rather than global offsets for the points
1454: Output Parameter:
1455: . gsection - The `PetscSection` for the global field layout
1457: Level: intermediate
1459: Notes:
1460: On each MPI process `gsection` inherits the chart of the `s` on that process.
1462: 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`.
1463: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1465: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCreateGlobalSectionCensored()`
1466: @*/
1467: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool usePermutation, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1468: {
1469: PetscSection gs;
1470: const PetscInt *pind = NULL;
1471: PetscInt *recv = NULL, *neg = NULL;
1472: PetscInt pStart, pEnd, p, dof, cdof, off, foff, globalOff = 0, nroots, nlocal, maxleaf;
1473: PetscInt numFields, f, numComponents;
1475: PetscFunctionBegin;
1481: PetscAssertPointer(gsection, 6);
1482: PetscCheck(s->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
1483: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs));
1484: PetscCall(PetscSectionGetNumFields(s, &numFields));
1485: if (numFields > 0) PetscCall(PetscSectionSetNumFields(gs, numFields));
1486: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1487: PetscCall(PetscSectionSetChart(gs, pStart, pEnd));
1488: gs->includesConstraints = includeConstraints;
1489: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1490: nlocal = nroots; /* The local/leaf space matches global/root space */
1491: /* Must allocate for all points visible to SF, which may be more than this section */
1492: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1493: PetscCall(PetscSFGetLeafRange(sf, NULL, &maxleaf));
1494: PetscCheck(nroots >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %" PetscInt_FMT " < pEnd %" PetscInt_FMT, nroots, pEnd);
1495: PetscCheck(maxleaf < nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %" PetscInt_FMT " >= nroots %" PetscInt_FMT, maxleaf, nroots);
1496: PetscCall(PetscMalloc2(nroots, &neg, nlocal, &recv));
1497: PetscCall(PetscArrayzero(neg, nroots));
1498: }
1499: /* Mark all local points with negative dof */
1500: for (p = pStart; p < pEnd; ++p) {
1501: PetscCall(PetscSectionGetDof(s, p, &dof));
1502: PetscCall(PetscSectionSetDof(gs, p, dof));
1503: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1504: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(gs, p, cdof));
1505: if (neg) neg[p] = -(dof + 1);
1506: }
1507: PetscCall(PetscSectionSetUpBC(gs));
1508: 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]));
1509: if (nroots >= 0) {
1510: PetscCall(PetscArrayzero(recv, nlocal));
1511: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1512: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1513: for (p = pStart; p < pEnd; ++p) {
1514: if (recv[p] < 0) {
1515: gs->atlasDof[p - pStart] = recv[p];
1516: PetscCall(PetscSectionGetDof(s, p, &dof));
1517: 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);
1518: }
1519: }
1520: }
1521: /* Calculate new sizes, get process offset, and calculate point offsets */
1522: if (usePermutation && s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1523: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1524: const PetscInt q = pind ? pind[p] : p;
1526: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1527: gs->atlasOff[q] = off;
1528: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1529: }
1530: if (!localOffsets) {
1531: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1532: globalOff -= off;
1533: }
1534: for (p = pStart, off = 0; p < pEnd; ++p) {
1535: gs->atlasOff[p - pStart] += globalOff;
1536: if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1537: }
1538: if (usePermutation && s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1539: /* Put in negative offsets for ghost points */
1540: if (nroots >= 0) {
1541: PetscCall(PetscArrayzero(recv, nlocal));
1542: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1543: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1544: for (p = pStart; p < pEnd; ++p) {
1545: if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1546: }
1547: }
1548: PetscCall(PetscFree2(neg, recv));
1549: /* Set field dofs/offsets/constraints */
1550: for (f = 0; f < numFields; ++f) {
1551: const char *name;
1553: gs->field[f]->includesConstraints = includeConstraints;
1554: PetscCall(PetscSectionGetFieldComponents(s, f, &numComponents));
1555: PetscCall(PetscSectionSetFieldComponents(gs, f, numComponents));
1556: PetscCall(PetscSectionGetFieldName(s, f, &name));
1557: PetscCall(PetscSectionSetFieldName(gs, f, name));
1558: }
1559: for (p = pStart; p < pEnd; ++p) {
1560: PetscCall(PetscSectionGetOffset(gs, p, &off));
1561: for (f = 0, foff = off; f < numFields; ++f) {
1562: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
1563: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetFieldConstraintDof(gs, p, f, cdof));
1564: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
1565: PetscCall(PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof));
1566: PetscCall(PetscSectionSetFieldOffset(gs, p, f, foff));
1567: PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof));
1568: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1569: }
1570: }
1571: for (f = 0; f < numFields; ++f) {
1572: PetscSection gfs = gs->field[f];
1574: PetscCall(PetscSectionSetUpBC(gfs));
1575: 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]));
1576: }
1577: gs->setup = PETSC_TRUE;
1578: PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view"));
1579: *gsection = gs;
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: /*@
1584: PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the globallayout using
1585: a local (sequential) `PetscSection` on each MPI process and an `PetscSF` describing the section point overlap.
1587: Input Parameters:
1588: + s - The `PetscSection` for the local field layout
1589: . sf - The `PetscSF` describing parallel layout of the section points
1590: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global vector will not possess constrained dofs
1591: . numExcludes - The number of exclusion ranges, this must have the same value on all MPI processes
1592: - 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
1594: Output Parameter:
1595: . gsection - The `PetscSection` for the global field layout
1597: Level: advanced
1599: Notes:
1600: On each MPI process `gsection` inherits the chart of the `s` on that process.
1602: 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`.
1603: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1605: This routine augments `PetscSectionCreateGlobalSection()` by allowing one to exclude certain ranges in the chart of the `PetscSection`
1607: Developer Notes:
1608: This is a terrible function name
1610: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1611: @*/
1612: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1613: {
1614: const PetscInt *pind = NULL;
1615: PetscInt *neg = NULL, *tmpOff = NULL;
1616: PetscInt pStart, pEnd, p, e, dof, cdof, off, globalOff = 0, nroots;
1618: PetscFunctionBegin;
1621: PetscAssertPointer(gsection, 6);
1622: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
1623: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1624: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
1625: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1626: if (nroots >= 0) {
1627: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
1628: PetscCall(PetscCalloc1(nroots, &neg));
1629: if (nroots > pEnd - pStart) {
1630: PetscCall(PetscCalloc1(nroots, &tmpOff));
1631: } else {
1632: tmpOff = &(*gsection)->atlasDof[-pStart];
1633: }
1634: }
1635: /* Mark ghost points with negative dof */
1636: for (p = pStart; p < pEnd; ++p) {
1637: for (e = 0; e < numExcludes; ++e) {
1638: if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1639: PetscCall(PetscSectionSetDof(*gsection, p, 0));
1640: break;
1641: }
1642: }
1643: if (e < numExcludes) continue;
1644: PetscCall(PetscSectionGetDof(s, p, &dof));
1645: PetscCall(PetscSectionSetDof(*gsection, p, dof));
1646: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1647: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
1648: if (neg) neg[p] = -(dof + 1);
1649: }
1650: PetscCall(PetscSectionSetUpBC(*gsection));
1651: if (nroots >= 0) {
1652: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1653: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1654: if (nroots > pEnd - pStart) {
1655: for (p = pStart; p < pEnd; ++p) {
1656: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1657: }
1658: }
1659: }
1660: /* Calculate new sizes, get process offset, and calculate point offsets */
1661: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1662: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1663: const PetscInt q = pind ? pind[p] : p;
1665: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1666: (*gsection)->atlasOff[q] = off;
1667: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1668: }
1669: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1670: globalOff -= off;
1671: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1672: (*gsection)->atlasOff[p] += globalOff;
1673: if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1674: }
1675: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1676: /* Put in negative offsets for ghost points */
1677: if (nroots >= 0) {
1678: if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1679: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1680: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1681: if (nroots > pEnd - pStart) {
1682: for (p = pStart; p < pEnd; ++p) {
1683: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1684: }
1685: }
1686: }
1687: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
1688: PetscCall(PetscFree(neg));
1689: PetscFunctionReturn(PETSC_SUCCESS);
1690: }
1692: /*@
1693: PetscSectionGetPointLayout - Get a `PetscLayout` for the points with nonzero dof counts of the unnamed default field within this `PetscSection`s local chart
1695: Collective
1697: Input Parameters:
1698: + comm - The `MPI_Comm`
1699: - s - The `PetscSection`
1701: Output Parameter:
1702: . layout - The point layout for the data that defines the section
1704: Level: advanced
1706: Notes:
1707: `PetscSectionGetValueLayout()` provides similar information but counting the total number of degrees of freedom on the MPI process (excluding constrained
1708: degrees of freedom).
1710: This count includes constrained degrees of freedom
1712: This is usually called on the default global section.
1714: Example:
1715: .vb
1716: The chart is [2,5), point 2 has 2 dof, point 3 has 0 dof, point 4 has 1 dof
1717: The local size of the `PetscLayout` is 2 since 2 points have a non-zero number of dof
1718: .ve
1720: Developer Notes:
1721: I find the names of these two functions extremely non-informative
1723: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1724: @*/
1725: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1726: {
1727: PetscInt pStart, pEnd, p, localSize = 0;
1729: PetscFunctionBegin;
1730: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1731: for (p = pStart; p < pEnd; ++p) {
1732: PetscInt dof;
1734: PetscCall(PetscSectionGetDof(s, p, &dof));
1735: if (dof >= 0) ++localSize;
1736: }
1737: PetscCall(PetscLayoutCreate(comm, layout));
1738: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1739: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1740: PetscCall(PetscLayoutSetUp(*layout));
1741: PetscFunctionReturn(PETSC_SUCCESS);
1742: }
1744: /*@
1745: PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs of a `PetscSection`
1747: Collective
1749: Input Parameters:
1750: + comm - The `MPI_Comm`
1751: - s - The `PetscSection`
1753: Output Parameter:
1754: . layout - The dof layout for the section
1756: Level: advanced
1758: Notes:
1759: `PetscSectionGetPointLayout()` provides similar information but only counting the number of points with nonzero degrees of freedom and
1760: including the constrained degrees of freedom
1762: This is usually called for the default global section.
1764: Example:
1765: .vb
1766: The chart is [2,5), point 2 has 4 dof (2 constrained), point 3 has 0 dof, point 4 has 1 dof (not constrained)
1767: The local size of the `PetscLayout` is 3 since there are 3 unconstrained degrees of freedom on this MPI process
1768: .ve
1770: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1771: @*/
1772: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1773: {
1774: PetscInt pStart, pEnd, p, localSize = 0;
1776: PetscFunctionBegin;
1778: PetscAssertPointer(layout, 3);
1779: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1780: for (p = pStart; p < pEnd; ++p) {
1781: PetscInt dof, cdof;
1783: PetscCall(PetscSectionGetDof(s, p, &dof));
1784: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1785: if (dof - cdof > 0) localSize += dof - cdof;
1786: }
1787: PetscCall(PetscLayoutCreate(comm, layout));
1788: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1789: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1790: PetscCall(PetscLayoutSetUp(*layout));
1791: PetscFunctionReturn(PETSC_SUCCESS);
1792: }
1794: /*@
1795: PetscSectionGetOffset - Return the offset into an array or `Vec` for the dof associated with the given point.
1797: Not Collective
1799: Input Parameters:
1800: + s - the `PetscSection`
1801: - point - the point
1803: Output Parameter:
1804: . offset - the offset
1806: Level: intermediate
1808: Notes:
1809: In a global section, `offset` will be negative for points not owned by this process.
1811: This is for the unnamed default field in the `PetscSection` not the named fields
1813: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1815: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetPointMajor()`
1816: @*/
1817: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1818: {
1819: PetscFunctionBegin;
1821: PetscAssertPointer(offset, 3);
1822: 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);
1823: *offset = s->atlasOff[point - s->pStart];
1824: PetscFunctionReturn(PETSC_SUCCESS);
1825: }
1827: /*@
1828: PetscSectionSetOffset - Set the offset into an array or `Vec` for the dof associated with the given point.
1830: Not Collective
1832: Input Parameters:
1833: + s - the `PetscSection`
1834: . point - the point
1835: - offset - the offset, these values may be negative indicating the values are off process
1837: Level: developer
1839: Note:
1840: The user usually does not call this function, but uses `PetscSectionSetUp()`
1842: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1843: @*/
1844: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1845: {
1846: PetscFunctionBegin;
1848: 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);
1849: s->atlasOff[point - s->pStart] = offset;
1850: PetscFunctionReturn(PETSC_SUCCESS);
1851: }
1853: /*@
1854: PetscSectionGetFieldOffset - Return the offset into an array or `Vec` for the field dof associated with the given point.
1856: Not Collective
1858: Input Parameters:
1859: + s - the `PetscSection`
1860: . point - the point
1861: - field - the field
1863: Output Parameter:
1864: . offset - the offset
1866: Level: intermediate
1868: Notes:
1869: In a global section, `offset` will be negative for points not owned by this process.
1871: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1873: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldPointOffset()`
1874: @*/
1875: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1876: {
1877: PetscFunctionBegin;
1879: PetscAssertPointer(offset, 4);
1880: PetscSectionCheckValidField(field, s->numFields);
1881: PetscCall(PetscSectionGetOffset(s->field[field], point, offset));
1882: PetscFunctionReturn(PETSC_SUCCESS);
1883: }
1885: /*@
1886: PetscSectionSetFieldOffset - Set the offset into an array or `Vec` for the dof associated with the given field at a point.
1888: Not Collective
1890: Input Parameters:
1891: + s - the `PetscSection`
1892: . point - the point
1893: . field - the field
1894: - offset - the offset, these values may be negative indicating the values are off process
1896: Level: developer
1898: Note:
1899: The user usually does not call this function, but uses `PetscSectionSetUp()`
1901: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1902: @*/
1903: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1904: {
1905: PetscFunctionBegin;
1907: PetscSectionCheckValidField(field, s->numFields);
1908: PetscCall(PetscSectionSetOffset(s->field[field], point, offset));
1909: PetscFunctionReturn(PETSC_SUCCESS);
1910: }
1912: /*@
1913: PetscSectionGetFieldPointOffset - Return the offset for the first field dof associated with the given point relative to the offset for that point for the
1914: unnamed default field's first dof
1916: Not Collective
1918: Input Parameters:
1919: + s - the `PetscSection`
1920: . point - the point
1921: - field - the field
1923: Output Parameter:
1924: . offset - the offset
1926: Level: advanced
1928: Note:
1929: This ignores constraints
1931: Example:
1932: .vb
1933: if PetscSectionSetPointMajor(s,PETSC_TRUE)
1934: The unnamed default field has 3 dof at `point`
1935: Field 0 has 2 dof at `point`
1936: Then PetscSectionGetFieldPointOffset(s,point,1,&offset) returns and offset of 5
1937: .ve
1939: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldOffset()`
1940: @*/
1941: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1942: {
1943: PetscInt off, foff;
1945: PetscFunctionBegin;
1947: PetscAssertPointer(offset, 4);
1948: PetscSectionCheckValidField(field, s->numFields);
1949: PetscCall(PetscSectionGetOffset(s, point, &off));
1950: PetscCall(PetscSectionGetOffset(s->field[field], point, &foff));
1951: *offset = foff - off;
1952: PetscFunctionReturn(PETSC_SUCCESS);
1953: }
1955: /*@
1956: PetscSectionGetOffsetRange - Return the full range of offsets [`start`, `end`) for a `PetscSection`
1958: Not Collective
1960: Input Parameter:
1961: . s - the `PetscSection`
1963: Output Parameters:
1964: + start - the minimum offset
1965: - end - one more than the maximum offset
1967: Level: intermediate
1969: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1970: @*/
1971: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1972: {
1973: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1975: PetscFunctionBegin;
1977: if (s->atlasOff) {
1978: os = s->atlasOff[0];
1979: oe = s->atlasOff[0];
1980: }
1981: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1982: for (p = 0; p < pEnd - pStart; ++p) {
1983: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1985: if (off >= 0) {
1986: os = PetscMin(os, off);
1987: oe = PetscMax(oe, off + dof);
1988: }
1989: }
1990: if (start) *start = os;
1991: if (end) *end = oe;
1992: PetscFunctionReturn(PETSC_SUCCESS);
1993: }
1995: /*@
1996: PetscSectionCreateSubsection - Create a new, smaller `PetscSection` composed of only selected fields
1998: Collective
2000: Input Parameters:
2001: + s - the `PetscSection`
2002: . len - the number of subfields
2003: - fields - the subfield numbers
2005: Output Parameter:
2006: . subs - the subsection
2008: Level: advanced
2010: Notes:
2011: The chart of `subs` is the same as the chart of `s`
2013: This will error if a fieldnumber is out of range
2015: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2016: @*/
2017: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
2018: {
2019: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
2021: PetscFunctionBegin;
2022: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2024: PetscAssertPointer(fields, 3);
2025: PetscAssertPointer(subs, 4);
2026: PetscCall(PetscSectionGetNumFields(s, &nF));
2027: 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);
2028: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2029: PetscCall(PetscSectionSetNumFields(*subs, len));
2030: for (f = 0; f < len; ++f) {
2031: const char *name = NULL;
2032: PetscInt numComp = 0;
2033: PetscSectionSym sym;
2035: PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
2036: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2037: PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
2038: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2039: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
2040: PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
2041: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2042: }
2043: PetscCall(PetscSectionGetFieldSym(s, fields[f], &sym));
2044: PetscCall(PetscSectionSetFieldSym(*subs, f, sym));
2045: }
2046: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2047: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2048: for (p = pStart; p < pEnd; ++p) {
2049: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
2051: for (f = 0; f < len; ++f) {
2052: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2053: PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
2054: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2055: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
2056: dof += fdof;
2057: cdof += cfdof;
2058: }
2059: PetscCall(PetscSectionSetDof(*subs, p, dof));
2060: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
2061: maxCdof = PetscMax(cdof, maxCdof);
2062: }
2063: PetscCall(PetscSectionSetUp(*subs));
2064: if (maxCdof) {
2065: PetscInt *indices;
2067: PetscCall(PetscMalloc1(maxCdof, &indices));
2068: for (p = pStart; p < pEnd; ++p) {
2069: PetscInt cdof;
2071: PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
2072: if (cdof) {
2073: const PetscInt *oldIndices = NULL;
2074: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
2076: for (f = 0; f < len; ++f) {
2077: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2078: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2079: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
2080: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
2081: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
2082: numConst += cfdof;
2083: fOff += fdof;
2084: }
2085: PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
2086: }
2087: }
2088: PetscCall(PetscFree(indices));
2089: }
2090: PetscFunctionReturn(PETSC_SUCCESS);
2091: }
2093: /*@
2094: PetscSectionCreateComponentSubsection - Create a new, smaller `PetscSection` composed of only selected components
2096: Collective
2098: Input Parameters:
2099: + s - the `PetscSection`
2100: . len - the number of components
2101: - comps - the component numbers
2103: Output Parameter:
2104: . subs - the subsection
2106: Level: advanced
2108: Notes:
2109: The chart of `subs` is the same as the chart of `s`
2111: This will error if the section has more than one field, or if a component number is out of range
2113: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2114: @*/
2115: PetscErrorCode PetscSectionCreateComponentSubsection(PetscSection s, PetscInt len, const PetscInt comps[], PetscSection *subs)
2116: {
2117: PetscSectionSym sym;
2118: const char *name = NULL;
2119: PetscInt Nf, pStart, pEnd;
2121: PetscFunctionBegin;
2122: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2124: PetscAssertPointer(comps, 3);
2125: PetscAssertPointer(subs, 4);
2126: PetscCall(PetscSectionGetNumFields(s, &Nf));
2127: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "This method can only handle one field, not %" PetscInt_FMT, Nf);
2128: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2129: PetscCall(PetscSectionSetNumFields(*subs, 1));
2130: PetscCall(PetscSectionGetFieldName(s, 0, &name));
2131: PetscCall(PetscSectionSetFieldName(*subs, 0, name));
2132: PetscCall(PetscSectionSetFieldComponents(*subs, 0, len));
2133: PetscCall(PetscSectionGetFieldSym(s, 0, &sym));
2134: PetscCall(PetscSectionSetFieldSym(*subs, 0, sym));
2135: for (PetscInt c = 0; c < len; ++c) {
2136: PetscCall(PetscSectionGetComponentName(s, 0, comps[c], &name));
2137: PetscCall(PetscSectionSetComponentName(*subs, 0, c, name));
2138: }
2139: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2140: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2141: for (PetscInt p = pStart; p < pEnd; ++p) {
2142: PetscInt dof, cdof, cfdof;
2144: PetscCall(PetscSectionGetDof(s, p, &dof));
2145: if (!dof) continue;
2146: PetscCall(PetscSectionGetFieldConstraintDof(s, p, 0, &cfdof));
2147: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2148: PetscCheck(!cdof && !cfdof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Component selection does not work with constraints");
2149: PetscCall(PetscSectionSetFieldDof(*subs, p, 0, len));
2150: PetscCall(PetscSectionSetDof(*subs, p, len));
2151: }
2152: PetscCall(PetscSectionSetUp(*subs));
2153: PetscFunctionReturn(PETSC_SUCCESS);
2154: }
2156: /*@
2157: PetscSectionCreateSupersection - Create a new, larger section composed of multiple `PetscSection`s
2159: Collective
2161: Input Parameters:
2162: + s - the input sections
2163: - len - the number of input sections
2165: Output Parameter:
2166: . supers - the supersection
2168: Level: advanced
2170: Notes:
2171: The section offsets now refer to a new, larger vector.
2173: Developer Notes:
2174: Needs to explain how the sections are composed
2176: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2177: @*/
2178: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2179: {
2180: PetscInt Nf = 0, f, pStart = PETSC_MAX_INT, pEnd = 0, p, maxCdof = 0, i;
2182: PetscFunctionBegin;
2183: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2184: for (i = 0; i < len; ++i) {
2185: PetscInt nf, pStarti, pEndi;
2187: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2188: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2189: pStart = PetscMin(pStart, pStarti);
2190: pEnd = PetscMax(pEnd, pEndi);
2191: Nf += nf;
2192: }
2193: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2194: PetscCall(PetscSectionSetNumFields(*supers, Nf));
2195: for (i = 0, f = 0; i < len; ++i) {
2196: PetscInt nf, fi, ci;
2198: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2199: for (fi = 0; fi < nf; ++fi, ++f) {
2200: const char *name = NULL;
2201: PetscInt numComp = 0;
2203: PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2204: PetscCall(PetscSectionSetFieldName(*supers, f, name));
2205: PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2206: PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2207: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2208: PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2209: PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2210: }
2211: }
2212: }
2213: PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2214: for (p = pStart; p < pEnd; ++p) {
2215: PetscInt dof = 0, cdof = 0;
2217: for (i = 0, f = 0; i < len; ++i) {
2218: PetscInt nf, fi, pStarti, pEndi;
2219: PetscInt fdof = 0, cfdof = 0;
2221: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2222: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2223: if ((p < pStarti) || (p >= pEndi)) continue;
2224: for (fi = 0; fi < nf; ++fi, ++f) {
2225: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2226: PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2227: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2228: if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2229: dof += fdof;
2230: cdof += cfdof;
2231: }
2232: }
2233: PetscCall(PetscSectionSetDof(*supers, p, dof));
2234: if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2235: maxCdof = PetscMax(cdof, maxCdof);
2236: }
2237: PetscCall(PetscSectionSetUp(*supers));
2238: if (maxCdof) {
2239: PetscInt *indices;
2241: PetscCall(PetscMalloc1(maxCdof, &indices));
2242: for (p = pStart; p < pEnd; ++p) {
2243: PetscInt cdof;
2245: PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2246: if (cdof) {
2247: PetscInt dof, numConst = 0, fOff = 0;
2249: for (i = 0, f = 0; i < len; ++i) {
2250: const PetscInt *oldIndices = NULL;
2251: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
2253: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2254: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2255: if ((p < pStarti) || (p >= pEndi)) continue;
2256: for (fi = 0; fi < nf; ++fi, ++f) {
2257: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2258: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2259: PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2260: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2261: PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2262: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2263: numConst += cfdof;
2264: }
2265: PetscCall(PetscSectionGetDof(s[i], p, &dof));
2266: fOff += dof;
2267: }
2268: PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2269: }
2270: }
2271: PetscCall(PetscFree(indices));
2272: }
2273: PetscFunctionReturn(PETSC_SUCCESS);
2274: }
2276: static PetscErrorCode PetscSectionCreateSubplexSection_Private(PetscSection s, IS subpointMap, PetscBool renumberPoints, PetscSection *subs)
2277: {
2278: const PetscInt *points = NULL, *indices = NULL;
2279: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
2281: PetscFunctionBegin;
2284: PetscAssertPointer(subs, 4);
2285: PetscCall(PetscSectionGetNumFields(s, &numFields));
2286: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2287: if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2288: for (f = 0; f < numFields; ++f) {
2289: const char *name = NULL;
2290: PetscInt numComp = 0;
2292: PetscCall(PetscSectionGetFieldName(s, f, &name));
2293: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2294: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2295: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2296: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2297: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2298: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2299: }
2300: }
2301: /* For right now, we do not try to squeeze the subchart */
2302: if (subpointMap) {
2303: PetscCall(ISGetSize(subpointMap, &numSubpoints));
2304: PetscCall(ISGetIndices(subpointMap, &points));
2305: }
2306: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2307: if (renumberPoints) {
2308: spStart = 0;
2309: spEnd = numSubpoints;
2310: } else {
2311: PetscCall(ISGetMinMax(subpointMap, &spStart, &spEnd));
2312: ++spEnd;
2313: }
2314: PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2315: for (p = pStart; p < pEnd; ++p) {
2316: PetscInt dof, cdof, fdof = 0, cfdof = 0;
2318: PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2319: if (subp < 0) continue;
2320: if (!renumberPoints) subp = p;
2321: for (f = 0; f < numFields; ++f) {
2322: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2323: PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2324: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2325: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2326: }
2327: PetscCall(PetscSectionGetDof(s, p, &dof));
2328: PetscCall(PetscSectionSetDof(*subs, subp, dof));
2329: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2330: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2331: }
2332: PetscCall(PetscSectionSetUp(*subs));
2333: /* Change offsets to original offsets */
2334: for (p = pStart; p < pEnd; ++p) {
2335: PetscInt off, foff = 0;
2337: PetscCall(PetscFindInt(p, numSubpoints, points, &subp));
2338: if (subp < 0) continue;
2339: if (!renumberPoints) subp = p;
2340: for (f = 0; f < numFields; ++f) {
2341: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2342: PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2343: }
2344: PetscCall(PetscSectionGetOffset(s, p, &off));
2345: PetscCall(PetscSectionSetOffset(*subs, subp, off));
2346: }
2347: /* Copy constraint indices */
2348: for (subp = spStart; subp < spEnd; ++subp) {
2349: PetscInt cdof;
2351: PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2352: if (cdof) {
2353: for (f = 0; f < numFields; ++f) {
2354: PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2355: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2356: }
2357: PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2358: PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2359: }
2360: }
2361: if (subpointMap) PetscCall(ISRestoreIndices(subpointMap, &points));
2362: PetscFunctionReturn(PETSC_SUCCESS);
2363: }
2365: /*@
2366: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2368: Collective
2370: Input Parameters:
2371: + s - the `PetscSection`
2372: - subpointMap - a sorted list of points in the original mesh which are in the submesh
2374: Output Parameter:
2375: . subs - the subsection
2377: Level: advanced
2379: Notes:
2380: 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))`
2382: Compare this with `PetscSectionCreateSubdomainSection()` that does not map the points numbers to start at zero but leaves them as before
2384: Developer Notes:
2385: 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`
2387: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2388: @*/
2389: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointMap, PetscSection *subs)
2390: {
2391: PetscFunctionBegin;
2392: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_TRUE, subs));
2393: PetscFunctionReturn(PETSC_SUCCESS);
2394: }
2396: /*@
2397: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2399: Collective
2401: Input Parameters:
2402: + s - the `PetscSection`
2403: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2405: Output Parameter:
2406: . subs - the subsection
2408: Level: advanced
2410: Notes:
2411: 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`
2412: is `[min(subpointMap),max(subpointMap)+1)`
2414: Compare this with `PetscSectionCreateSubmeshSection()` that maps the point numbers to start at zero
2416: Developer Notes:
2417: 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`
2419: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2420: @*/
2421: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2422: {
2423: PetscFunctionBegin;
2424: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_FALSE, subs));
2425: PetscFunctionReturn(PETSC_SUCCESS);
2426: }
2428: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2429: {
2430: PetscInt p;
2431: PetscMPIInt rank;
2433: PetscFunctionBegin;
2434: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2435: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2436: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2437: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2438: if (s->bc && s->bc->atlasDof[p] > 0) {
2439: PetscInt b;
2440: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2441: if (s->bcIndices) {
2442: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2443: }
2444: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2445: } else {
2446: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dim %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2447: }
2448: }
2449: PetscCall(PetscViewerFlush(viewer));
2450: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2451: if (s->sym) {
2452: PetscCall(PetscViewerASCIIPushTab(viewer));
2453: PetscCall(PetscSectionSymView(s->sym, viewer));
2454: PetscCall(PetscViewerASCIIPopTab(viewer));
2455: }
2456: PetscFunctionReturn(PETSC_SUCCESS);
2457: }
2459: /*@C
2460: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2462: Collective
2464: Input Parameters:
2465: + A - the `PetscSection` object to view
2466: . obj - Optional object that provides the options prefix used for the options
2467: - name - command line option
2469: Level: intermediate
2471: Note:
2472: See `PetscObjectViewFromOptions()` for available values of `PetscViewer` and `PetscViewerFormat`
2474: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2475: @*/
2476: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2477: {
2478: PetscFunctionBegin;
2480: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2481: PetscFunctionReturn(PETSC_SUCCESS);
2482: }
2484: /*@C
2485: PetscSectionView - Views a `PetscSection`
2487: Collective
2489: Input Parameters:
2490: + s - the `PetscSection` object to view
2491: - viewer - the viewer
2493: Level: beginner
2495: Note:
2496: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2497: distribution independent data, such as dofs, offsets, constraint dofs,
2498: and constraint indices. Points that have negative dofs, for instance,
2499: are not saved as they represent points owned by other processes.
2500: Point numbering and rank assignment is currently not stored.
2501: The saved section can be loaded with `PetscSectionLoad()`.
2503: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2504: @*/
2505: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2506: {
2507: PetscBool isascii, ishdf5;
2508: PetscInt f;
2510: PetscFunctionBegin;
2512: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2514: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2515: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2516: if (isascii) {
2517: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2518: if (s->numFields) {
2519: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2520: for (f = 0; f < s->numFields; ++f) {
2521: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " \"%s\" with %" PetscInt_FMT " components\n", f, s->fieldNames[f], s->numFieldComponents[f]));
2522: PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2523: }
2524: } else {
2525: PetscCall(PetscSectionView_ASCII(s, viewer));
2526: }
2527: } else if (ishdf5) {
2528: #if PetscDefined(HAVE_HDF5)
2529: PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2530: #else
2531: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2532: #endif
2533: }
2534: PetscFunctionReturn(PETSC_SUCCESS);
2535: }
2537: /*@C
2538: PetscSectionLoad - Loads a `PetscSection`
2540: Collective
2542: Input Parameters:
2543: + s - the `PetscSection` object to load
2544: - viewer - the viewer
2546: Level: beginner
2548: Note:
2549: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2550: a section saved with `PetscSectionView()`. The number of processes
2551: used here (N) does not need to be the same as that used when saving.
2552: After calling this function, the chart of s on rank i will be set
2553: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2554: saved section points.
2556: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2557: @*/
2558: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2559: {
2560: PetscBool ishdf5;
2562: PetscFunctionBegin;
2565: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2566: if (ishdf5) {
2567: #if PetscDefined(HAVE_HDF5)
2568: PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2569: PetscFunctionReturn(PETSC_SUCCESS);
2570: #else
2571: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2572: #endif
2573: } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2574: }
2576: /*@
2577: PetscSectionResetClosurePermutation - Remove any existing closure permutation
2579: Input Parameter:
2580: . section - The `PetscSection`
2582: Level: intermediate
2584: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionSetClosureIndex()`, `PetscSectionReset()`
2585: @*/
2586: PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2587: {
2588: PetscSectionClosurePermVal clVal;
2590: PetscFunctionBegin;
2591: if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2592: kh_foreach_value(section->clHash, clVal, {
2593: PetscCall(PetscFree(clVal.perm));
2594: PetscCall(PetscFree(clVal.invPerm));
2595: });
2596: kh_destroy(ClPerm, section->clHash);
2597: section->clHash = NULL;
2598: PetscFunctionReturn(PETSC_SUCCESS);
2599: }
2601: /*@
2602: PetscSectionReset - Frees all section data, the section is then as if `PetscSectionCreate()` had just been called.
2604: Not Collective
2606: Input Parameter:
2607: . s - the `PetscSection`
2609: Level: beginner
2611: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2612: @*/
2613: PetscErrorCode PetscSectionReset(PetscSection s)
2614: {
2615: PetscInt f, c;
2617: PetscFunctionBegin;
2619: for (f = 0; f < s->numFields; ++f) {
2620: PetscCall(PetscSectionDestroy(&s->field[f]));
2621: PetscCall(PetscFree(s->fieldNames[f]));
2622: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2623: PetscCall(PetscFree(s->compNames[f]));
2624: }
2625: PetscCall(PetscFree(s->numFieldComponents));
2626: PetscCall(PetscFree(s->fieldNames));
2627: PetscCall(PetscFree(s->compNames));
2628: PetscCall(PetscFree(s->field));
2629: PetscCall(PetscSectionDestroy(&s->bc));
2630: PetscCall(PetscFree(s->bcIndices));
2631: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2632: PetscCall(PetscSectionDestroy(&s->clSection));
2633: PetscCall(ISDestroy(&s->clPoints));
2634: PetscCall(ISDestroy(&s->perm));
2635: PetscCall(PetscBTDestroy(&s->blockStarts));
2636: PetscCall(PetscSectionResetClosurePermutation(s));
2637: PetscCall(PetscSectionSymDestroy(&s->sym));
2638: PetscCall(PetscSectionDestroy(&s->clSection));
2639: PetscCall(ISDestroy(&s->clPoints));
2640: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2641: s->pStart = -1;
2642: s->pEnd = -1;
2643: s->maxDof = 0;
2644: s->setup = PETSC_FALSE;
2645: s->numFields = 0;
2646: s->clObj = NULL;
2647: PetscFunctionReturn(PETSC_SUCCESS);
2648: }
2650: /*@
2651: PetscSectionDestroy - Frees a `PetscSection`
2653: Not Collective
2655: Input Parameter:
2656: . s - the `PetscSection`
2658: Level: beginner
2660: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2661: @*/
2662: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2663: {
2664: PetscFunctionBegin;
2665: if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2667: if (--((PetscObject)*s)->refct > 0) {
2668: *s = NULL;
2669: PetscFunctionReturn(PETSC_SUCCESS);
2670: }
2671: PetscCall(PetscSectionReset(*s));
2672: PetscCall(PetscHeaderDestroy(s));
2673: PetscFunctionReturn(PETSC_SUCCESS);
2674: }
2676: static PetscErrorCode VecIntGetValuesSection_Private(const PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2677: {
2678: const PetscInt p = point - s->pStart;
2680: PetscFunctionBegin;
2682: *values = &baseArray[s->atlasOff[p]];
2683: PetscFunctionReturn(PETSC_SUCCESS);
2684: }
2686: static PetscErrorCode VecIntSetValuesSection_Private(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2687: {
2688: PetscInt *array;
2689: const PetscInt p = point - s->pStart;
2690: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2691: PetscInt cDim = 0;
2693: PetscFunctionBegin;
2695: PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2696: array = &baseArray[s->atlasOff[p]];
2697: if (!cDim) {
2698: if (orientation >= 0) {
2699: const PetscInt dim = s->atlasDof[p];
2700: PetscInt i;
2702: if (mode == INSERT_VALUES) {
2703: for (i = 0; i < dim; ++i) array[i] = values ? values[i] : i;
2704: } else {
2705: for (i = 0; i < dim; ++i) array[i] += values[i];
2706: }
2707: } else {
2708: PetscInt offset = 0;
2709: PetscInt j = -1, field, i;
2711: for (field = 0; field < s->numFields; ++field) {
2712: const PetscInt dim = s->field[field]->atlasDof[p];
2714: for (i = dim - 1; i >= 0; --i) array[++j] = values ? values[i + offset] : i + offset;
2715: offset += dim;
2716: }
2717: }
2718: } else {
2719: if (orientation >= 0) {
2720: const PetscInt dim = s->atlasDof[p];
2721: PetscInt cInd = 0, i;
2722: const PetscInt *cDof;
2724: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2725: if (mode == INSERT_VALUES) {
2726: for (i = 0; i < dim; ++i) {
2727: if ((cInd < cDim) && (i == cDof[cInd])) {
2728: ++cInd;
2729: continue;
2730: }
2731: array[i] = values ? values[i] : i;
2732: }
2733: } else {
2734: for (i = 0; i < dim; ++i) {
2735: if ((cInd < cDim) && (i == cDof[cInd])) {
2736: ++cInd;
2737: continue;
2738: }
2739: array[i] += values[i];
2740: }
2741: }
2742: } else {
2743: const PetscInt *cDof;
2744: PetscInt offset = 0;
2745: PetscInt cOffset = 0;
2746: PetscInt j = 0, field;
2748: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2749: for (field = 0; field < s->numFields; ++field) {
2750: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2751: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2752: const PetscInt sDim = dim - tDim;
2753: PetscInt cInd = 0, i, k;
2755: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2756: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2757: ++cInd;
2758: continue;
2759: }
2760: array[j] = values ? values[k] : k;
2761: }
2762: offset += dim;
2763: cOffset += dim - tDim;
2764: }
2765: }
2766: }
2767: PetscFunctionReturn(PETSC_SUCCESS);
2768: }
2770: /*@C
2771: PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs
2773: Not Collective
2775: Input Parameter:
2776: . s - The `PetscSection`
2778: Output Parameter:
2779: . hasConstraints - flag indicating that the section has constrained dofs
2781: Level: intermediate
2783: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2784: @*/
2785: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2786: {
2787: PetscFunctionBegin;
2789: PetscAssertPointer(hasConstraints, 2);
2790: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2791: PetscFunctionReturn(PETSC_SUCCESS);
2792: }
2794: /*@C
2795: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained for a given point
2797: Not Collective
2799: Input Parameters:
2800: + s - The `PetscSection`
2801: - point - The point
2803: Output Parameter:
2804: . indices - The constrained dofs
2806: Level: intermediate
2808: Fortran Notes:
2809: Use `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`
2811: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2812: @*/
2813: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt **indices)
2814: {
2815: PetscFunctionBegin;
2817: if (s->bc) {
2818: PetscCall(VecIntGetValuesSection_Private(s->bcIndices, s->bc, point, indices));
2819: } else *indices = NULL;
2820: PetscFunctionReturn(PETSC_SUCCESS);
2821: }
2823: /*@C
2824: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2826: Not Collective
2828: Input Parameters:
2829: + s - The `PetscSection`
2830: . point - The point
2831: - indices - The constrained dofs
2833: Level: intermediate
2835: Fortran Notes:
2836: Use `PetscSectionSetConstraintIndicesF90()`
2838: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2839: @*/
2840: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2841: {
2842: PetscFunctionBegin;
2844: if (s->bc) {
2845: const PetscInt dof = s->atlasDof[point];
2846: const PetscInt cdof = s->bc->atlasDof[point];
2847: PetscInt d;
2849: if (indices)
2850: for (d = 0; d < cdof; ++d) 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]);
2851: PetscCall(VecIntSetValuesSection_Private(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
2852: }
2853: PetscFunctionReturn(PETSC_SUCCESS);
2854: }
2856: /*@C
2857: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2859: Not Collective
2861: Input Parameters:
2862: + s - The `PetscSection`
2863: . field - The field number
2864: - point - The point
2866: Output Parameter:
2867: . indices - The constrained dofs sorted in ascending order
2869: Level: intermediate
2871: Note:
2872: The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by `PetscSectionGetConstraintDof()`.
2874: Fortran Notes:
2875: Use `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`
2877: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2878: @*/
2879: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2880: {
2881: PetscFunctionBegin;
2883: PetscAssertPointer(indices, 4);
2884: PetscSectionCheckValidField(field, s->numFields);
2885: PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
2886: PetscFunctionReturn(PETSC_SUCCESS);
2887: }
2889: /*@C
2890: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2892: Not Collective
2894: Input Parameters:
2895: + s - The `PetscSection`
2896: . point - The point
2897: . field - The field number
2898: - indices - The constrained dofs
2900: Level: intermediate
2902: Fortran Notes:
2903: Use `PetscSectionSetFieldConstraintIndicesF90()`
2905: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2906: @*/
2907: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2908: {
2909: PetscFunctionBegin;
2911: PetscSectionCheckValidField(field, s->numFields);
2912: PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
2913: PetscFunctionReturn(PETSC_SUCCESS);
2914: }
2916: /*@
2917: PetscSectionPermute - Reorder the section according to the input point permutation
2919: Collective
2921: Input Parameters:
2922: + section - The `PetscSection` object
2923: - permutation - The point permutation, old point p becomes new point perm[p]
2925: Output Parameter:
2926: . sectionNew - The permuted `PetscSection`
2928: Level: intermediate
2930: Note:
2931: The data and the access to the data via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` are both changed in `sectionNew`
2933: Compare to `PetscSectionSetPermutation()`
2935: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
2936: @*/
2937: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2938: {
2939: PetscSection s = section, sNew;
2940: const PetscInt *perm;
2941: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2943: PetscFunctionBegin;
2946: PetscAssertPointer(sectionNew, 3);
2947: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
2948: PetscCall(PetscSectionGetNumFields(s, &numFields));
2949: if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
2950: for (f = 0; f < numFields; ++f) {
2951: const char *name;
2952: PetscInt numComp;
2954: PetscCall(PetscSectionGetFieldName(s, f, &name));
2955: PetscCall(PetscSectionSetFieldName(sNew, f, name));
2956: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2957: PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
2958: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2959: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2960: PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
2961: }
2962: }
2963: PetscCall(ISGetLocalSize(permutation, &numPoints));
2964: PetscCall(ISGetIndices(permutation, &perm));
2965: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2966: PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
2967: PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
2968: for (p = pStart; p < pEnd; ++p) {
2969: PetscInt dof, cdof;
2971: PetscCall(PetscSectionGetDof(s, p, &dof));
2972: PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
2973: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2974: if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
2975: for (f = 0; f < numFields; ++f) {
2976: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
2977: PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
2978: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2979: if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
2980: }
2981: }
2982: PetscCall(PetscSectionSetUp(sNew));
2983: for (p = pStart; p < pEnd; ++p) {
2984: const PetscInt *cind;
2985: PetscInt cdof;
2987: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2988: if (cdof) {
2989: PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
2990: PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
2991: }
2992: for (f = 0; f < numFields; ++f) {
2993: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2994: if (cdof) {
2995: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
2996: PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
2997: }
2998: }
2999: }
3000: PetscCall(ISRestoreIndices(permutation, &perm));
3001: *sectionNew = sNew;
3002: PetscFunctionReturn(PETSC_SUCCESS);
3003: }
3005: /*@
3006: PetscSectionSetClosureIndex - Create an internal data structure to speed up closure queries.
3008: Collective
3010: Input Parameters:
3011: + section - The `PetscSection`
3012: . obj - A `PetscObject` which serves as the key for this index
3013: . clSection - `PetscSection` giving the size of the closure of each point
3014: - clPoints - `IS` giving the points in each closure
3016: Level: advanced
3018: Note:
3019: This function creates an internal map from each point to its closure. We compress out closure points with no dofs in this section.
3021: Developer Notes:
3022: The information provided here is completely opaque
3024: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
3025: @*/
3026: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
3027: {
3028: PetscFunctionBegin;
3032: if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
3033: section->clObj = obj;
3034: PetscCall(PetscObjectReference((PetscObject)clSection));
3035: PetscCall(PetscObjectReference((PetscObject)clPoints));
3036: PetscCall(PetscSectionDestroy(§ion->clSection));
3037: PetscCall(ISDestroy(§ion->clPoints));
3038: section->clSection = clSection;
3039: section->clPoints = clPoints;
3040: PetscFunctionReturn(PETSC_SUCCESS);
3041: }
3043: /*@
3044: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
3046: Collective
3048: Input Parameters:
3049: + section - The `PetscSection`
3050: - obj - A `PetscObject` which serves as the key for this index
3052: Output Parameters:
3053: + clSection - `PetscSection` giving the size of the closure of each point
3054: - clPoints - `IS` giving the points in each closure
3056: Level: advanced
3058: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3059: @*/
3060: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
3061: {
3062: PetscFunctionBegin;
3063: if (section->clObj == obj) {
3064: if (clSection) *clSection = section->clSection;
3065: if (clPoints) *clPoints = section->clPoints;
3066: } else {
3067: if (clSection) *clSection = NULL;
3068: if (clPoints) *clPoints = NULL;
3069: }
3070: PetscFunctionReturn(PETSC_SUCCESS);
3071: }
3073: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
3074: {
3075: PetscInt i;
3076: khiter_t iter;
3077: int new_entry;
3078: PetscSectionClosurePermKey key = {depth, clSize};
3079: PetscSectionClosurePermVal *val;
3081: PetscFunctionBegin;
3082: if (section->clObj != obj) {
3083: PetscCall(PetscSectionDestroy(§ion->clSection));
3084: PetscCall(ISDestroy(§ion->clPoints));
3085: }
3086: section->clObj = obj;
3087: if (!section->clHash) PetscCall(PetscClPermCreate(§ion->clHash));
3088: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
3089: val = &kh_val(section->clHash, iter);
3090: if (!new_entry) {
3091: PetscCall(PetscFree(val->perm));
3092: PetscCall(PetscFree(val->invPerm));
3093: }
3094: if (mode == PETSC_COPY_VALUES) {
3095: PetscCall(PetscMalloc1(clSize, &val->perm));
3096: PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
3097: } else if (mode == PETSC_OWN_POINTER) {
3098: val->perm = clPerm;
3099: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
3100: PetscCall(PetscMalloc1(clSize, &val->invPerm));
3101: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
3102: PetscFunctionReturn(PETSC_SUCCESS);
3103: }
3105: /*@
3106: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3108: Not Collective
3110: Input Parameters:
3111: + section - The `PetscSection`
3112: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3113: . depth - Depth of points on which to apply the given permutation
3114: - perm - Permutation of the cell dof closure
3116: Level: intermediate
3118: Notes:
3119: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
3120: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
3121: each topology and degree.
3123: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
3124: exotic/enriched spaces on mixed topology meshes.
3126: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
3127: @*/
3128: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
3129: {
3130: const PetscInt *clPerm = NULL;
3131: PetscInt clSize = 0;
3133: PetscFunctionBegin;
3134: if (perm) {
3135: PetscCall(ISGetLocalSize(perm, &clSize));
3136: PetscCall(ISGetIndices(perm, &clPerm));
3137: }
3138: PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
3139: if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
3140: PetscFunctionReturn(PETSC_SUCCESS);
3141: }
3143: static PetscErrorCode PetscSectionGetClosurePermutation_Private(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3144: {
3145: PetscFunctionBegin;
3146: if (section->clObj == obj) {
3147: PetscSectionClosurePermKey k = {depth, size};
3148: PetscSectionClosurePermVal v;
3150: PetscCall(PetscClPermGet(section->clHash, k, &v));
3151: if (perm) *perm = v.perm;
3152: } else {
3153: if (perm) *perm = NULL;
3154: }
3155: PetscFunctionReturn(PETSC_SUCCESS);
3156: }
3158: /*@
3159: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3161: Not Collective
3163: Input Parameters:
3164: + section - The `PetscSection`
3165: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
3166: . depth - Depth stratum on which to obtain closure permutation
3167: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3169: Output Parameter:
3170: . perm - The dof closure permutation
3172: Level: intermediate
3174: Note:
3175: The user must destroy the `IS` that is returned.
3177: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3178: @*/
3179: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3180: {
3181: const PetscInt *clPerm = NULL;
3183: PetscFunctionBegin;
3184: PetscCall(PetscSectionGetClosurePermutation_Private(section, obj, depth, clSize, &clPerm));
3185: PetscCheck(clPerm, PetscObjectComm((PetscObject)obj), PETSC_ERR_ARG_WRONG, "There is no closure permutation associated with this object for depth %" PetscInt_FMT " of size %" PetscInt_FMT, depth, clSize);
3186: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3187: PetscFunctionReturn(PETSC_SUCCESS);
3188: }
3190: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3191: {
3192: PetscFunctionBegin;
3193: if (section->clObj == obj && section->clHash) {
3194: PetscSectionClosurePermKey k = {depth, size};
3195: PetscSectionClosurePermVal v;
3196: PetscCall(PetscClPermGet(section->clHash, k, &v));
3197: if (perm) *perm = v.invPerm;
3198: } else {
3199: if (perm) *perm = NULL;
3200: }
3201: PetscFunctionReturn(PETSC_SUCCESS);
3202: }
3204: /*@
3205: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
3207: Not Collective
3209: Input Parameters:
3210: + section - The `PetscSection`
3211: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3212: . depth - Depth stratum on which to obtain closure permutation
3213: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3215: Output Parameter:
3216: . perm - The dof closure permutation
3218: Level: intermediate
3220: Note:
3221: The user must destroy the `IS` that is returned.
3223: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3224: @*/
3225: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3226: {
3227: const PetscInt *clPerm = NULL;
3229: PetscFunctionBegin;
3230: PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3231: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3232: PetscFunctionReturn(PETSC_SUCCESS);
3233: }
3235: /*@
3236: PetscSectionGetField - Get the `PetscSection` associated with a single field
3238: Input Parameters:
3239: + s - The `PetscSection`
3240: - field - The field number
3242: Output Parameter:
3243: . subs - The `PetscSection` for the given field, note the chart of `subs` is not set
3245: Level: intermediate
3247: Note:
3248: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
3250: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3251: @*/
3252: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3253: {
3254: PetscFunctionBegin;
3256: PetscAssertPointer(subs, 3);
3257: PetscSectionCheckValidField(field, s->numFields);
3258: *subs = s->field[field];
3259: PetscFunctionReturn(PETSC_SUCCESS);
3260: }
3262: PetscClassId PETSC_SECTION_SYM_CLASSID;
3263: PetscFunctionList PetscSectionSymList = NULL;
3265: /*@
3266: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
3268: Collective
3270: Input Parameter:
3271: . comm - the MPI communicator
3273: Output Parameter:
3274: . sym - pointer to the new set of symmetries
3276: Level: developer
3278: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3279: @*/
3280: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3281: {
3282: PetscFunctionBegin;
3283: PetscAssertPointer(sym, 2);
3284: PetscCall(ISInitializePackage());
3285: PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3286: PetscFunctionReturn(PETSC_SUCCESS);
3287: }
3289: /*@C
3290: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
3292: Collective
3294: Input Parameters:
3295: + sym - The section symmetry object
3296: - method - The name of the section symmetry type
3298: Level: developer
3300: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3301: @*/
3302: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3303: {
3304: PetscErrorCode (*r)(PetscSectionSym);
3305: PetscBool match;
3307: PetscFunctionBegin;
3309: PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3310: if (match) PetscFunctionReturn(PETSC_SUCCESS);
3312: PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3313: PetscCheck(r, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3314: PetscTryTypeMethod(sym, destroy);
3315: sym->ops->destroy = NULL;
3317: PetscCall((*r)(sym));
3318: PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3319: PetscFunctionReturn(PETSC_SUCCESS);
3320: }
3322: /*@C
3323: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
3325: Not Collective
3327: Input Parameter:
3328: . sym - The section symmetry
3330: Output Parameter:
3331: . type - The index set type name
3333: Level: developer
3335: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3336: @*/
3337: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3338: {
3339: PetscFunctionBegin;
3341: PetscAssertPointer(type, 2);
3342: *type = ((PetscObject)sym)->type_name;
3343: PetscFunctionReturn(PETSC_SUCCESS);
3344: }
3346: /*@C
3347: PetscSectionSymRegister - Registers a new section symmetry implementation
3349: Not Collective
3351: Input Parameters:
3352: + sname - The name of a new user-defined creation routine
3353: - function - The creation routine itself
3355: Level: developer
3357: Notes:
3358: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
3360: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3361: @*/
3362: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3363: {
3364: PetscFunctionBegin;
3365: PetscCall(ISInitializePackage());
3366: PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3367: PetscFunctionReturn(PETSC_SUCCESS);
3368: }
3370: /*@
3371: PetscSectionSymDestroy - Destroys a section symmetry.
3373: Collective
3375: Input Parameter:
3376: . sym - the section symmetry
3378: Level: developer
3380: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`
3381: @*/
3382: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3383: {
3384: SymWorkLink link, next;
3386: PetscFunctionBegin;
3387: if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3389: if (--((PetscObject)*sym)->refct > 0) {
3390: *sym = NULL;
3391: PetscFunctionReturn(PETSC_SUCCESS);
3392: }
3393: PetscTryTypeMethod(*sym, destroy);
3394: PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3395: for (link = (*sym)->workin; link; link = next) {
3396: PetscInt **perms = (PetscInt **)link->perms;
3397: PetscScalar **rots = (PetscScalar **)link->rots;
3398: PetscCall(PetscFree2(perms, rots));
3399: next = link->next;
3400: PetscCall(PetscFree(link));
3401: }
3402: (*sym)->workin = NULL;
3403: PetscCall(PetscHeaderDestroy(sym));
3404: PetscFunctionReturn(PETSC_SUCCESS);
3405: }
3407: /*@C
3408: PetscSectionSymView - Displays a section symmetry
3410: Collective
3412: Input Parameters:
3413: + sym - the index set
3414: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3416: Level: developer
3418: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3419: @*/
3420: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3421: {
3422: PetscFunctionBegin;
3424: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3426: PetscCheckSameComm(sym, 1, viewer, 2);
3427: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3428: PetscTryTypeMethod(sym, view, viewer);
3429: PetscFunctionReturn(PETSC_SUCCESS);
3430: }
3432: /*@
3433: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3435: Collective
3437: Input Parameters:
3438: + section - the section describing data layout
3439: - sym - the symmetry describing the affect of orientation on the access of the data
3441: Level: developer
3443: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3444: @*/
3445: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3446: {
3447: PetscFunctionBegin;
3449: PetscCall(PetscSectionSymDestroy(§ion->sym));
3450: if (sym) {
3452: PetscCheckSameComm(section, 1, sym, 2);
3453: PetscCall(PetscObjectReference((PetscObject)sym));
3454: }
3455: section->sym = sym;
3456: PetscFunctionReturn(PETSC_SUCCESS);
3457: }
3459: /*@
3460: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3462: Not Collective
3464: Input Parameter:
3465: . section - the section describing data layout
3467: Output Parameter:
3468: . sym - the symmetry describing the affect of orientation on the access of the data, provided previously by `PetscSectionSetSym()`
3470: Level: developer
3472: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3473: @*/
3474: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3475: {
3476: PetscFunctionBegin;
3478: *sym = section->sym;
3479: PetscFunctionReturn(PETSC_SUCCESS);
3480: }
3482: /*@
3483: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3485: Collective
3487: Input Parameters:
3488: + section - the section describing data layout
3489: . field - the field number
3490: - sym - the symmetry describing the affect of orientation on the access of the data
3492: Level: developer
3494: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3495: @*/
3496: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3497: {
3498: PetscFunctionBegin;
3500: PetscSectionCheckValidField(field, section->numFields);
3501: PetscCall(PetscSectionSetSym(section->field[field], sym));
3502: PetscFunctionReturn(PETSC_SUCCESS);
3503: }
3505: /*@
3506: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3508: Collective
3510: Input Parameters:
3511: + section - the section describing data layout
3512: - field - the field number
3514: Output Parameter:
3515: . sym - the symmetry describing the affect of orientation on the access of the data
3517: Level: developer
3519: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3520: @*/
3521: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3522: {
3523: PetscFunctionBegin;
3525: PetscSectionCheckValidField(field, section->numFields);
3526: *sym = section->field[field]->sym;
3527: PetscFunctionReturn(PETSC_SUCCESS);
3528: }
3530: /*@C
3531: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3533: Not Collective
3535: Input Parameters:
3536: + section - the section
3537: . numPoints - the number of points
3538: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3539: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3540: context, see `DMPlexGetConeOrientation()`).
3542: Output Parameters:
3543: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3544: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3545: identity).
3547: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3548: .vb
3549: const PetscInt **perms;
3550: const PetscScalar **rots;
3551: PetscInt lOffset;
3553: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3554: for (i = 0, lOffset = 0; i < numPoints; i++) {
3555: PetscInt point = points[2*i], dof, sOffset;
3556: const PetscInt *perm = perms ? perms[i] : NULL;
3557: const PetscScalar *rot = rots ? rots[i] : NULL;
3559: PetscSectionGetDof(section,point,&dof);
3560: PetscSectionGetOffset(section,point,&sOffset);
3562: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3563: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3564: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3565: lOffset += dof;
3566: }
3567: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3568: .ve
3570: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3571: .vb
3572: const PetscInt **perms;
3573: const PetscScalar **rots;
3574: PetscInt lOffset;
3576: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3577: for (i = 0, lOffset = 0; i < numPoints; i++) {
3578: PetscInt point = points[2*i], dof, sOffset;
3579: const PetscInt *perm = perms ? perms[i] : NULL;
3580: const PetscScalar *rot = rots ? rots[i] : NULL;
3582: PetscSectionGetDof(section,point,&dof);
3583: PetscSectionGetOffset(section,point,&sOff);
3585: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3586: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3587: offset += dof;
3588: }
3589: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3590: .ve
3592: Level: developer
3594: Notes:
3595: `PetscSectionSetSym()` must have been previously called to provide the symmetries to the `PetscSection`
3597: Use `PetscSectionRestorePointSyms()` when finished with the data
3599: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3600: @*/
3601: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3602: {
3603: PetscSectionSym sym;
3605: PetscFunctionBegin;
3607: if (numPoints) PetscAssertPointer(points, 3);
3608: if (perms) *perms = NULL;
3609: if (rots) *rots = NULL;
3610: sym = section->sym;
3611: if (sym && (perms || rots)) {
3612: SymWorkLink link;
3614: if (sym->workin) {
3615: link = sym->workin;
3616: sym->workin = sym->workin->next;
3617: } else {
3618: PetscCall(PetscNew(&link));
3619: }
3620: if (numPoints > link->numPoints) {
3621: PetscInt **perms = (PetscInt **)link->perms;
3622: PetscScalar **rots = (PetscScalar **)link->rots;
3623: PetscCall(PetscFree2(perms, rots));
3624: PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3625: link->numPoints = numPoints;
3626: }
3627: link->next = sym->workout;
3628: sym->workout = link;
3629: PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3630: PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3631: PetscUseTypeMethod(sym, getpoints, section, numPoints, points, link->perms, link->rots);
3632: if (perms) *perms = link->perms;
3633: if (rots) *rots = link->rots;
3634: }
3635: PetscFunctionReturn(PETSC_SUCCESS);
3636: }
3638: /*@C
3639: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3641: Not Collective
3643: Input Parameters:
3644: + section - the section
3645: . numPoints - the number of points
3646: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3647: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3648: context, see `DMPlexGetConeOrientation()`).
3649: . perms - The permutations for the given orientations: set to `NULL` at conclusion
3650: - rots - The field rotations symmetries for the given orientations: set to `NULL` at conclusion
3652: Level: developer
3654: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3655: @*/
3656: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3657: {
3658: PetscSectionSym sym;
3660: PetscFunctionBegin;
3662: sym = section->sym;
3663: if (sym && (perms || rots)) {
3664: SymWorkLink *p, link;
3666: for (p = &sym->workout; (link = *p); p = &link->next) {
3667: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3668: *p = link->next;
3669: link->next = sym->workin;
3670: sym->workin = link;
3671: if (perms) *perms = NULL;
3672: if (rots) *rots = NULL;
3673: PetscFunctionReturn(PETSC_SUCCESS);
3674: }
3675: }
3676: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3677: }
3678: PetscFunctionReturn(PETSC_SUCCESS);
3679: }
3681: /*@C
3682: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3684: Not Collective
3686: Input Parameters:
3687: + section - the section
3688: . field - the field of the section
3689: . numPoints - the number of points
3690: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3691: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3692: context, see `DMPlexGetConeOrientation()`).
3694: Output Parameters:
3695: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3696: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3697: identity).
3699: Level: developer
3701: Notes:
3702: `PetscSectionSetFieldSym()` must have been previously called to provide the symmetries to the `PetscSection`
3704: Use `PetscSectionRestoreFieldPointSyms()` when finished with the data
3706: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3707: @*/
3708: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3709: {
3710: PetscFunctionBegin;
3712: 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);
3713: PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3714: PetscFunctionReturn(PETSC_SUCCESS);
3715: }
3717: /*@C
3718: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3720: Not Collective
3722: Input Parameters:
3723: + section - the section
3724: . field - the field number
3725: . numPoints - the number of points
3726: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3727: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3728: context, see `DMPlexGetConeOrientation()`).
3729: . perms - The permutations for the given orientations: set to NULL at conclusion
3730: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3732: Level: developer
3734: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3735: @*/
3736: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3737: {
3738: PetscFunctionBegin;
3740: 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);
3741: PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3742: PetscFunctionReturn(PETSC_SUCCESS);
3743: }
3745: /*@
3746: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3748: Not Collective
3750: Input Parameter:
3751: . sym - the `PetscSectionSym`
3753: Output Parameter:
3754: . nsym - the equivalent symmetries
3756: Level: developer
3758: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3759: @*/
3760: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3761: {
3762: PetscFunctionBegin;
3765: PetscTryTypeMethod(sym, copy, nsym);
3766: PetscFunctionReturn(PETSC_SUCCESS);
3767: }
3769: /*@
3770: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3772: Collective
3774: Input Parameters:
3775: + sym - the `PetscSectionSym`
3776: - migrationSF - the distribution map from roots to leaves
3778: Output Parameter:
3779: . dsym - the redistributed symmetries
3781: Level: developer
3783: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3784: @*/
3785: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3786: {
3787: PetscFunctionBegin;
3790: PetscAssertPointer(dsym, 3);
3791: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3792: PetscFunctionReturn(PETSC_SUCCESS);
3793: }
3795: /*@
3796: PetscSectionGetUseFieldOffsets - Get the flag indicating if field offsets are used directly in a global section, rather than just the point offset
3798: Not Collective
3800: Input Parameter:
3801: . s - the global `PetscSection`
3803: Output Parameter:
3804: . flg - the flag
3806: Level: developer
3808: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3809: @*/
3810: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3811: {
3812: PetscFunctionBegin;
3814: *flg = s->useFieldOff;
3815: PetscFunctionReturn(PETSC_SUCCESS);
3816: }
3818: /*@
3819: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3821: Not Collective
3823: Input Parameters:
3824: + s - the global `PetscSection`
3825: - flg - the flag
3827: Level: developer
3829: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3830: @*/
3831: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3832: {
3833: PetscFunctionBegin;
3835: s->useFieldOff = flg;
3836: PetscFunctionReturn(PETSC_SUCCESS);
3837: }
3839: #define PetscSectionExpandPoints_Loop(TYPE) \
3840: do { \
3841: PetscInt i, n, o0, o1, size; \
3842: TYPE *a0 = (TYPE *)origArray, *a1; \
3843: PetscCall(PetscSectionGetStorageSize(s, &size)); \
3844: PetscCall(PetscMalloc1(size, &a1)); \
3845: for (i = 0; i < npoints; i++) { \
3846: PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
3847: PetscCall(PetscSectionGetOffset(s, i, &o1)); \
3848: PetscCall(PetscSectionGetDof(s, i, &n)); \
3849: PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n *unitsize)); \
3850: } \
3851: *newArray = (void *)a1; \
3852: } while (0)
3854: /*@
3855: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3857: Not Collective
3859: Input Parameters:
3860: + origSection - the `PetscSection` describing the layout of the array
3861: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
3862: . origArray - the array; its size must be equal to the storage size of `origSection`
3863: - points - `IS` with points to extract; its indices must lie in the chart of `origSection`
3865: Output Parameters:
3866: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3867: - newArray - the array of the extracted DOFs; its size is the storage size of `newSection`
3869: Level: developer
3871: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3872: @*/
3873: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3874: {
3875: PetscSection s;
3876: const PetscInt *points_;
3877: PetscInt i, n, npoints, pStart, pEnd;
3878: PetscMPIInt unitsize;
3880: PetscFunctionBegin;
3882: PetscAssertPointer(origArray, 3);
3884: if (newSection) PetscAssertPointer(newSection, 5);
3885: if (newArray) PetscAssertPointer(newArray, 6);
3886: PetscCallMPI(MPI_Type_size(dataType, &unitsize));
3887: PetscCall(ISGetLocalSize(points, &npoints));
3888: PetscCall(ISGetIndices(points, &points_));
3889: PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
3890: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
3891: PetscCall(PetscSectionSetChart(s, 0, npoints));
3892: for (i = 0; i < npoints; i++) {
3893: 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);
3894: PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
3895: PetscCall(PetscSectionSetDof(s, i, n));
3896: }
3897: PetscCall(PetscSectionSetUp(s));
3898: if (newArray) {
3899: if (dataType == MPIU_INT) {
3900: PetscSectionExpandPoints_Loop(PetscInt);
3901: } else if (dataType == MPIU_SCALAR) {
3902: PetscSectionExpandPoints_Loop(PetscScalar);
3903: } else if (dataType == MPIU_REAL) {
3904: PetscSectionExpandPoints_Loop(PetscReal);
3905: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3906: }
3907: if (newSection) {
3908: *newSection = s;
3909: } else {
3910: PetscCall(PetscSectionDestroy(&s));
3911: }
3912: PetscCall(ISRestoreIndices(points, &points_));
3913: PetscFunctionReturn(PETSC_SUCCESS);
3914: }