Actual source code: section.c
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsf.h>
8: PetscClassId PETSC_SECTION_CLASSID;
10: /*@
11: PetscSectionCreate - Allocates a `PetscSection` and sets the map contents to the default.
13: Collective
15: Input Parameters:
16: + comm - the MPI communicator
17: - s - pointer to the section
19: Level: beginner
21: Notes:
22: Typical calling sequence
23: .vb
24: PetscSectionCreate(MPI_Comm,PetscSection *);!
25: PetscSectionSetNumFields(PetscSection, numFields);
26: PetscSectionSetChart(PetscSection,low,high);
27: PetscSectionSetDof(PetscSection,point,numdof);
28: PetscSectionSetUp(PetscSection);
29: PetscSectionGetOffset(PetscSection,point,PetscInt *);
30: PetscSectionDestroy(PetscSection);
31: .ve
33: The `PetscSection` object and methods are intended to be used in the PETSc `Vec` and `Mat` implementations. The indices returned by the `PetscSection` are appropriate for the kind of `Vec` it is associated with. For example, if the vector being indexed is a local vector, we call the section a local section. If the section indexes a global vector, we call it a global section. For parallel vectors, like global vectors, we use negative indices to indicate dofs owned by other processes.
35: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionDestroy()`, `PetscSectionCreateGlobalSection()`
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
39: PetscFunctionBegin;
40: PetscAssertPointer(s, 2);
41: PetscCall(ISInitializePackage());
43: PetscCall(PetscHeaderCreate(*s, PETSC_SECTION_CLASSID, "PetscSection", "Section", "IS", comm, PetscSectionDestroy, PetscSectionView));
44: (*s)->pStart = -1;
45: (*s)->pEnd = -1;
46: (*s)->perm = NULL;
47: (*s)->pointMajor = PETSC_TRUE;
48: (*s)->includesConstraints = PETSC_TRUE;
49: (*s)->atlasDof = NULL;
50: (*s)->atlasOff = NULL;
51: (*s)->bc = NULL;
52: (*s)->bcIndices = NULL;
53: (*s)->setup = PETSC_FALSE;
54: (*s)->numFields = 0;
55: (*s)->fieldNames = NULL;
56: (*s)->field = NULL;
57: (*s)->useFieldOff = PETSC_FALSE;
58: (*s)->compNames = NULL;
59: (*s)->clObj = NULL;
60: (*s)->clHash = NULL;
61: (*s)->clSection = NULL;
62: (*s)->clPoints = NULL;
63: PetscCall(PetscSectionInvalidateMaxDof_Internal(*s));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*@
68: PetscSectionCopy - Creates a shallow (if possible) copy of the `PetscSection`
70: Collective
72: Input Parameter:
73: . section - the `PetscSection`
75: Output Parameter:
76: . newSection - the copy
78: Level: intermediate
80: Developer Notes:
81: What exactly does shallow mean in this context?
83: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
84: @*/
85: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
86: {
87: PetscFunctionBegin;
90: PetscCall(PetscSectionCopy_Internal(section, newSection, NULL));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: PetscErrorCode PetscSectionCopy_Internal(PetscSection section, PetscSection newSection, PetscBT constrained_dofs)
95: {
96: PetscSectionSym sym;
97: IS perm;
98: PetscInt numFields, f, c, pStart, pEnd, p;
100: PetscFunctionBegin;
103: PetscCall(PetscSectionReset(newSection));
104: PetscCall(PetscSectionGetNumFields(section, &numFields));
105: if (numFields) PetscCall(PetscSectionSetNumFields(newSection, numFields));
106: for (f = 0; f < numFields; ++f) {
107: const char *fieldName = NULL, *compName = NULL;
108: PetscInt numComp = 0;
110: PetscCall(PetscSectionGetFieldName(section, f, &fieldName));
111: PetscCall(PetscSectionSetFieldName(newSection, f, fieldName));
112: PetscCall(PetscSectionGetFieldComponents(section, f, &numComp));
113: PetscCall(PetscSectionSetFieldComponents(newSection, f, numComp));
114: for (c = 0; c < numComp; ++c) {
115: PetscCall(PetscSectionGetComponentName(section, f, c, &compName));
116: PetscCall(PetscSectionSetComponentName(newSection, f, c, compName));
117: }
118: PetscCall(PetscSectionGetFieldSym(section, f, &sym));
119: PetscCall(PetscSectionSetFieldSym(newSection, f, sym));
120: }
121: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
122: PetscCall(PetscSectionSetChart(newSection, pStart, pEnd));
123: PetscCall(PetscSectionGetPermutation(section, &perm));
124: PetscCall(PetscSectionSetPermutation(newSection, perm));
125: PetscCall(PetscSectionGetSym(section, &sym));
126: PetscCall(PetscSectionSetSym(newSection, sym));
127: for (p = pStart; p < pEnd; ++p) {
128: PetscInt dof, cdof, fcdof = 0;
129: PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));
131: PetscCall(PetscSectionGetDof(section, p, &dof));
132: PetscCall(PetscSectionSetDof(newSection, p, dof));
133: if (force_constrained) cdof = dof;
134: else PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
135: if (cdof) PetscCall(PetscSectionSetConstraintDof(newSection, p, cdof));
136: for (f = 0; f < numFields; ++f) {
137: PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
138: PetscCall(PetscSectionSetFieldDof(newSection, p, f, dof));
139: if (cdof) {
140: if (force_constrained) fcdof = dof;
141: else PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
142: if (fcdof) PetscCall(PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof));
143: }
144: }
145: }
146: PetscCall(PetscSectionSetUp(newSection));
147: for (p = pStart; p < pEnd; ++p) {
148: PetscInt off, cdof, fcdof = 0;
149: const PetscInt *cInd;
150: PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));
152: /* Must set offsets in case they do not agree with the prefix sums */
153: PetscCall(PetscSectionGetOffset(section, p, &off));
154: PetscCall(PetscSectionSetOffset(newSection, p, off));
155: PetscCall(PetscSectionGetConstraintDof(newSection, p, &cdof));
156: if (cdof) {
157: if (force_constrained) cInd = NULL;
158: else PetscCall(PetscSectionGetConstraintIndices(section, p, &cInd));
159: PetscCall(PetscSectionSetConstraintIndices(newSection, p, cInd));
160: for (f = 0; f < numFields; ++f) {
161: PetscCall(PetscSectionGetFieldOffset(section, p, f, &off));
162: PetscCall(PetscSectionSetFieldOffset(newSection, p, f, off));
163: PetscCall(PetscSectionGetFieldConstraintDof(newSection, p, f, &fcdof));
164: if (fcdof) {
165: if (force_constrained) cInd = NULL;
166: else PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &cInd));
167: PetscCall(PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd));
168: }
169: }
170: }
171: }
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: /*@
176: PetscSectionClone - Creates a shallow (if possible) copy of the `PetscSection`
178: Collective
180: Input Parameter:
181: . section - the `PetscSection`
183: Output Parameter:
184: . newSection - the copy
186: Level: beginner
188: Developer Notes:
189: With standard PETSc terminology this should be called `PetscSectionDuplicate()`
191: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
192: @*/
193: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
194: {
195: PetscFunctionBegin;
197: PetscAssertPointer(newSection, 2);
198: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection));
199: PetscCall(PetscSectionCopy(section, *newSection));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: PetscSectionSetFromOptions - sets parameters in a `PetscSection` from the options database
206: Collective
208: Input Parameter:
209: . s - the `PetscSection`
211: Options Database Key:
212: . -petscsection_point_major - `PETSC_TRUE` for point-major order
214: Level: intermediate
216: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
217: @*/
218: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
219: {
220: PetscFunctionBegin;
222: PetscObjectOptionsBegin((PetscObject)s);
223: PetscCall(PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL));
224: /* process any options handlers added with PetscObjectAddOptionsHandler() */
225: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject));
226: PetscOptionsEnd();
227: PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view"));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: PetscSectionCompare - Compares two sections
234: Collective
236: Input Parameters:
237: + s1 - the first `PetscSection`
238: - s2 - the second `PetscSection`
240: Output Parameter:
241: . congruent - `PETSC_TRUE` if the two sections are congruent, `PETSC_FALSE` otherwise
243: Level: intermediate
245: Note:
246: Field names are disregarded.
248: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
249: @*/
250: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
251: {
252: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
253: const PetscInt *idx1, *idx2;
254: IS perm1, perm2;
255: PetscBool flg;
256: PetscMPIInt mflg;
258: PetscFunctionBegin;
261: PetscAssertPointer(congruent, 3);
262: flg = PETSC_FALSE;
264: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg));
265: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
266: *congruent = PETSC_FALSE;
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: PetscCall(PetscSectionGetChart(s1, &pStart, &pEnd));
271: PetscCall(PetscSectionGetChart(s2, &n1, &n2));
272: if (pStart != n1 || pEnd != n2) goto not_congruent;
274: PetscCall(PetscSectionGetPermutation(s1, &perm1));
275: PetscCall(PetscSectionGetPermutation(s2, &perm2));
276: if (perm1 && perm2) {
277: PetscCall(ISEqual(perm1, perm2, congruent));
278: if (!*congruent) goto not_congruent;
279: } else if (perm1 != perm2) goto not_congruent;
281: for (p = pStart; p < pEnd; ++p) {
282: PetscCall(PetscSectionGetOffset(s1, p, &n1));
283: PetscCall(PetscSectionGetOffset(s2, p, &n2));
284: if (n1 != n2) goto not_congruent;
286: PetscCall(PetscSectionGetDof(s1, p, &n1));
287: PetscCall(PetscSectionGetDof(s2, p, &n2));
288: if (n1 != n2) goto not_congruent;
290: PetscCall(PetscSectionGetConstraintDof(s1, p, &ncdof));
291: PetscCall(PetscSectionGetConstraintDof(s2, p, &n2));
292: if (ncdof != n2) goto not_congruent;
294: PetscCall(PetscSectionGetConstraintIndices(s1, p, &idx1));
295: PetscCall(PetscSectionGetConstraintIndices(s2, p, &idx2));
296: PetscCall(PetscArraycmp(idx1, idx2, ncdof, congruent));
297: if (!*congruent) goto not_congruent;
298: }
300: PetscCall(PetscSectionGetNumFields(s1, &nfields));
301: PetscCall(PetscSectionGetNumFields(s2, &n2));
302: if (nfields != n2) goto not_congruent;
304: for (f = 0; f < nfields; ++f) {
305: PetscCall(PetscSectionGetFieldComponents(s1, f, &n1));
306: PetscCall(PetscSectionGetFieldComponents(s2, f, &n2));
307: if (n1 != n2) goto not_congruent;
309: for (p = pStart; p < pEnd; ++p) {
310: PetscCall(PetscSectionGetFieldOffset(s1, p, f, &n1));
311: PetscCall(PetscSectionGetFieldOffset(s2, p, f, &n2));
312: if (n1 != n2) goto not_congruent;
314: PetscCall(PetscSectionGetFieldDof(s1, p, f, &n1));
315: PetscCall(PetscSectionGetFieldDof(s2, p, f, &n2));
316: if (n1 != n2) goto not_congruent;
318: PetscCall(PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof));
319: PetscCall(PetscSectionGetFieldConstraintDof(s2, p, f, &n2));
320: if (nfcdof != n2) goto not_congruent;
322: PetscCall(PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1));
323: PetscCall(PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2));
324: PetscCall(PetscArraycmp(idx1, idx2, nfcdof, congruent));
325: if (!*congruent) goto not_congruent;
326: }
327: }
329: flg = PETSC_TRUE;
330: not_congruent:
331: PetscCallMPI(MPIU_Allreduce(&flg, congruent, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1)));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: /*@
336: PetscSectionGetNumFields - Returns the number of fields in a `PetscSection`, or 0 if no fields were defined.
338: Not Collective
340: Input Parameter:
341: . s - the `PetscSection`
343: Output Parameter:
344: . numFields - the number of fields defined, or 0 if none were defined
346: Level: intermediate
348: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
349: @*/
350: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
351: {
352: PetscFunctionBegin;
354: PetscAssertPointer(numFields, 2);
355: *numFields = s->numFields;
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: /*@
360: PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`
362: Not Collective
364: Input Parameters:
365: + s - the `PetscSection`
366: - numFields - the number of fields
368: Level: intermediate
370: Notes:
371: Calling this destroys all the information in the `PetscSection` including the chart.
373: You must call `PetscSectionSetChart()` after calling this.
375: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetNumFields()`, `PetscSectionSetChart()`, `PetscSectionReset()`
376: @*/
377: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
378: {
379: PetscInt f;
381: PetscFunctionBegin;
383: PetscCheck(numFields > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %" PetscInt_FMT " must be positive", numFields);
384: PetscCall(PetscSectionReset(s));
386: s->numFields = numFields;
387: PetscCall(PetscMalloc1(s->numFields, &s->numFieldComponents));
388: PetscCall(PetscMalloc1(s->numFields, &s->fieldNames));
389: PetscCall(PetscMalloc1(s->numFields, &s->compNames));
390: PetscCall(PetscMalloc1(s->numFields, &s->field));
391: for (f = 0; f < s->numFields; ++f) {
392: char name[64];
394: s->numFieldComponents[f] = 1;
396: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]));
397: PetscCall(PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f));
398: PetscCall(PetscStrallocpy(name, &s->fieldNames[f]));
399: PetscCall(PetscSNPrintf(name, 64, "Component_0"));
400: PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]));
401: PetscCall(PetscStrallocpy(name, &s->compNames[f][0]));
402: }
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: /*@
407: PetscSectionGetFieldName - Returns the name of a field in the `PetscSection`
409: Not Collective
411: Input Parameters:
412: + s - the `PetscSection`
413: - field - the field number
415: Output Parameter:
416: . fieldName - the field name
418: Level: intermediate
420: Note:
421: Will error if the field number is out of range
423: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
424: @*/
425: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
426: {
427: PetscFunctionBegin;
429: PetscAssertPointer(fieldName, 3);
430: PetscSectionCheckValidField(field, s->numFields);
431: *fieldName = s->fieldNames[field];
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: /*@
436: PetscSectionSetFieldName - Sets the name of a field in the `PetscSection`
438: Not Collective
440: Input Parameters:
441: + s - the `PetscSection`
442: . field - the field number
443: - fieldName - the field name
445: Level: intermediate
447: Note:
448: Will error if the field number is out of range
450: .seealso: [PetscSection](ch_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
451: @*/
452: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
453: {
454: PetscFunctionBegin;
456: if (fieldName) PetscAssertPointer(fieldName, 3);
457: PetscSectionCheckValidField(field, s->numFields);
458: PetscCall(PetscFree(s->fieldNames[field]));
459: PetscCall(PetscStrallocpy(fieldName, &s->fieldNames[field]));
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: /*@
464: PetscSectionGetComponentName - Gets the name of a field component in the `PetscSection`
466: Not Collective
468: Input Parameters:
469: + s - the `PetscSection`
470: . field - the field number
471: - comp - the component number
473: Output Parameter:
474: . compName - the component name
476: Level: intermediate
478: Note:
479: Will error if the field or component number do not exist
481: Developer Notes:
482: The function name should have Field in it since they are field components.
484: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
485: `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
486: @*/
487: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
488: {
489: PetscFunctionBegin;
491: PetscAssertPointer(compName, 4);
492: PetscSectionCheckValidField(field, s->numFields);
493: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
494: *compName = s->compNames[field][comp];
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /*@
499: PetscSectionSetComponentName - Sets the name of a field component in the `PetscSection`
501: Not Collective
503: Input Parameters:
504: + s - the `PetscSection`
505: . field - the field number
506: . comp - the component number
507: - compName - the component name
509: Level: advanced
511: Note:
512: Will error if the field or component number do not exist
514: Developer Notes:
515: The function name should have Field in it since they are field components.
517: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
518: `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
519: @*/
520: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
521: {
522: PetscFunctionBegin;
524: if (compName) PetscAssertPointer(compName, 4);
525: PetscSectionCheckValidField(field, s->numFields);
526: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
527: PetscCall(PetscFree(s->compNames[field][comp]));
528: PetscCall(PetscStrallocpy(compName, &s->compNames[field][comp]));
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: /*@
533: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
535: Not Collective
537: Input Parameters:
538: + s - the `PetscSection`
539: - field - the field number
541: Output Parameter:
542: . numComp - the number of field components
544: Level: advanced
546: Developer Notes:
547: This function is misnamed. There is a Num in `PetscSectionGetNumFields()` but not in this name
549: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`,
550: `PetscSectionSetComponentName()`, `PetscSectionGetComponentName()`
551: @*/
552: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
553: {
554: PetscFunctionBegin;
556: PetscAssertPointer(numComp, 3);
557: PetscSectionCheckValidField(field, s->numFields);
558: *numComp = s->numFieldComponents[field];
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: /*@
563: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
565: Not Collective
567: Input Parameters:
568: + s - the `PetscSection`
569: . field - the field number
570: - numComp - the number of field components
572: Level: advanced
574: Note:
575: This number can be different than the values set with `PetscSectionSetFieldDof()`. It can be used to indicate the number of
576: components in the field of the underlying physical model which may be different than the number of degrees of freedom needed
577: 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
578: an face based model for velocity (where the velocity normal to the face is stored) there is only 1 dof for each face point.
580: The value set with this function are not needed or used in `PetscSectionSetUp()`.
582: Developer Notes:
583: This function is misnamed. There is a Num in `PetscSectionSetNumFields()` but not in this name
585: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionSetComponentName()`,
586: `PetscSectionGetComponentName()`, `PetscSectionGetNumFields()`
587: @*/
588: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
589: {
590: PetscInt c;
592: PetscFunctionBegin;
594: PetscSectionCheckValidField(field, s->numFields);
595: if (s->compNames) {
596: for (c = 0; c < s->numFieldComponents[field]; ++c) PetscCall(PetscFree(s->compNames[field][c]));
597: PetscCall(PetscFree(s->compNames[field]));
598: }
600: s->numFieldComponents[field] = numComp;
601: if (numComp) {
602: PetscCall(PetscMalloc1(numComp, &s->compNames[field]));
603: for (c = 0; c < numComp; ++c) {
604: char name[64];
606: PetscCall(PetscSNPrintf(name, 64, "%" PetscInt_FMT, c));
607: PetscCall(PetscStrallocpy(name, &s->compNames[field][c]));
608: }
609: }
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: /*@
614: PetscSectionGetChart - Returns the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
616: Not Collective
618: Input Parameter:
619: . s - the `PetscSection`
621: Output Parameters:
622: + pStart - the first point
623: - pEnd - one past the last point
625: Level: intermediate
627: Note:
628: The chart may be thought of as the bounds on the points (indices) one may use to index into numerical data that is associated with
629: the `PetscSection` data layout.
631: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
632: @*/
633: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
634: {
635: PetscFunctionBegin;
637: if (pStart) *pStart = s->pStart;
638: if (pEnd) *pEnd = s->pEnd;
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: /*@
643: PetscSectionSetChart - Sets the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
645: Not Collective
647: Input Parameters:
648: + s - the `PetscSection`
649: . pStart - the first `point`
650: - pEnd - one past the last point, `pStart` $ \le $ `pEnd`
652: Level: intermediate
654: Notes:
655: The chart may be thought of as the bounds on the points (indices) one may use to index into numerical data that is associated with
656: the `PetscSection` data layout.
658: The charts on different MPI processes may (and often do) overlap
660: If you intend to use `PetscSectionSetNumFields()` it must be called before this call.
662: The chart for all fields created with `PetscSectionSetNumFields()` is the same as this chart.
664: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`, `PetscSectionSetNumFields()`
665: @*/
666: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
667: {
668: PetscInt f;
670: PetscFunctionBegin;
672: PetscCheck(pEnd >= pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Chart pEnd %" PetscInt_FMT " cannot be smaller than chart pStart %" PetscInt_FMT, pEnd, pStart);
673: if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(PETSC_SUCCESS);
674: /* Cannot Reset() because it destroys field information */
675: s->setup = PETSC_FALSE;
676: PetscCall(PetscSectionDestroy(&s->bc));
677: PetscCall(PetscFree(s->bcIndices));
678: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
680: s->pStart = pStart;
681: s->pEnd = pEnd;
682: PetscCall(PetscMalloc2(pEnd - pStart, &s->atlasDof, pEnd - pStart, &s->atlasOff));
683: PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart));
684: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd));
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: /*@
689: PetscSectionGetPermutation - Returns the permutation of [0, `pEnd` - `pStart`) or `NULL` that was set with `PetscSectionSetPermutation()`
691: Not Collective
693: Input Parameter:
694: . s - the `PetscSection`
696: Output Parameter:
697: . perm - The permutation as an `IS`
699: Level: intermediate
701: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
702: @*/
703: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
704: {
705: PetscFunctionBegin;
707: if (perm) {
708: PetscAssertPointer(perm, 2);
709: *perm = s->perm;
710: }
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: /*@
715: PetscSectionSetPermutation - Sets a permutation of the chart for this section, [0, `pEnd` - `pStart`), which determines the order to store the `PetscSection` information
717: Not Collective
719: Input Parameters:
720: + s - the `PetscSection`
721: - perm - the permutation of points
723: Level: intermediate
725: Notes:
726: The permutation must be provided before `PetscSectionSetUp()`.
728: The data in the `PetscSection` are permuted but the access via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` is not changed
730: Compare to `PetscSectionPermute()`
732: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetUp()`, `PetscSectionGetPermutation()`, `PetscSectionPermute()`, `PetscSectionCreate()`
733: @*/
734: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
735: {
736: PetscFunctionBegin;
739: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
740: if (s->perm != perm) {
741: PetscCall(ISDestroy(&s->perm));
742: if (perm) {
743: s->perm = perm;
744: PetscCall(PetscObjectReference((PetscObject)s->perm));
745: }
746: }
747: PetscFunctionReturn(PETSC_SUCCESS);
748: }
750: /*@C
751: PetscSectionGetBlockStarts - Returns a table indicating which points start new blocks
753: Not Collective, No Fortran Support
755: Input Parameter:
756: . s - the `PetscSection`
758: Output Parameter:
759: . blockStarts - The `PetscBT` with a 1 for each point that begins a block
761: Notes:
762: The table is on [0, `pEnd` - `pStart`).
764: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
766: Level: intermediate
768: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
769: @*/
770: PetscErrorCode PetscSectionGetBlockStarts(PetscSection s, PetscBT *blockStarts)
771: {
772: PetscFunctionBegin;
774: if (blockStarts) {
775: PetscAssertPointer(blockStarts, 2);
776: *blockStarts = s->blockStarts;
777: }
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: /*@C
782: PetscSectionSetBlockStarts - Sets a table indicating which points start new blocks
784: Not Collective, No Fortran Support
786: Input Parameters:
787: + s - the `PetscSection`
788: - blockStarts - The `PetscBT` with a 1 for each point that begins a block
790: Level: intermediate
792: Notes:
793: 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.
795: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
797: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionGetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
798: @*/
799: PetscErrorCode PetscSectionSetBlockStarts(PetscSection s, PetscBT blockStarts)
800: {
801: PetscFunctionBegin;
803: if (s->blockStarts != blockStarts) {
804: PetscCall(PetscBTDestroy(&s->blockStarts));
805: s->blockStarts = blockStarts;
806: }
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: /*@
811: PetscSectionGetPointMajor - Returns the flag for dof ordering, `PETSC_TRUE` if it is point major, `PETSC_FALSE` if it is field major
813: Not Collective
815: Input Parameter:
816: . s - the `PetscSection`
818: Output Parameter:
819: . pm - the flag for point major ordering
821: Level: intermediate
823: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetPointMajor()`
824: @*/
825: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
826: {
827: PetscFunctionBegin;
829: PetscAssertPointer(pm, 2);
830: *pm = s->pointMajor;
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: /*@
835: PetscSectionSetPointMajor - Sets the flag for dof ordering, `PETSC_TRUE` for point major, otherwise it will be field major
837: Not Collective
839: Input Parameters:
840: + s - the `PetscSection`
841: - pm - the flag for point major ordering
843: Level: intermediate
845: Note:
846: 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.
848: Point major order means the degrees of freedom are stored as follows
849: .vb
850: all the degrees of freedom for each point are stored contiguously, one point after another (respecting a permutation set with `PetscSectionSetPermutation()`)
851: for each point
852: the degrees of freedom for each field (starting with the unnamed default field) are listed in order by field
853: .ve
855: Field major order means the degrees of freedom are stored as follows
856: .vb
857: all degrees of freedom for each field (including the unnamed default field) are stored contiguously, one field after another
858: for each field (started with unnamed default field)
859: the degrees of freedom for each point are listed in order by point (respecting a permutation set with `PetscSectionSetPermutation()`)
860: .ve
862: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`, `PetscSectionSetPermutation()`
863: @*/
864: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
865: {
866: PetscFunctionBegin;
868: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
869: s->pointMajor = pm;
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*@
874: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets in the `PetscSection`.
875: The value is set with `PetscSectionSetIncludesConstraints()`
877: Not Collective
879: Input Parameter:
880: . s - the `PetscSection`
882: Output Parameter:
883: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
885: Level: intermediate
887: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
888: @*/
889: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
890: {
891: PetscFunctionBegin;
893: PetscAssertPointer(includesConstraints, 2);
894: *includesConstraints = s->includesConstraints;
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: /*@
899: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
901: Not Collective
903: Input Parameters:
904: + s - the `PetscSection`
905: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
907: Level: intermediate
909: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
910: @*/
911: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
912: {
913: PetscFunctionBegin;
915: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
916: s->includesConstraints = includesConstraints;
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: /*@
921: PetscSectionGetDof - Return the total number of degrees of freedom associated with a given point.
923: Not Collective
925: Input Parameters:
926: + s - the `PetscSection`
927: - point - the point
929: Output Parameter:
930: . numDof - the number of dof
932: Level: intermediate
934: Notes:
935: In a global section, this size will be negative for points not owned by this process.
937: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
939: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
940: @*/
941: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
942: {
943: PetscFunctionBeginHot;
945: PetscAssertPointer(numDof, 3);
946: 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);
947: *numDof = s->atlasDof[point - s->pStart];
948: PetscFunctionReturn(PETSC_SUCCESS);
949: }
951: /*@
952: PetscSectionSetDof - Sets the total number of degrees of freedom associated with a given point.
954: Not Collective
956: Input Parameters:
957: + s - the `PetscSection`
958: . point - the point
959: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
961: Level: intermediate
963: Note:
964: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
966: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
967: @*/
968: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
969: {
970: PetscFunctionBegin;
972: 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);
973: s->atlasDof[point - s->pStart] = numDof;
974: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
975: PetscFunctionReturn(PETSC_SUCCESS);
976: }
978: /*@
979: PetscSectionAddDof - Adds to the total number of degrees of freedom associated with a given point.
981: Not Collective
983: Input Parameters:
984: + s - the `PetscSection`
985: . point - the point
986: - numDof - the number of additional dof
988: Level: intermediate
990: Note:
991: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
993: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
994: @*/
995: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
996: {
997: PetscFunctionBeginHot;
999: 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);
1000: PetscCheck(numDof >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numDof %" PetscInt_FMT " should not be negative", numDof);
1001: s->atlasDof[point - s->pStart] += numDof;
1002: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: /*@
1007: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
1009: Not Collective
1011: Input Parameters:
1012: + s - the `PetscSection`
1013: . point - the point
1014: - field - the field
1016: Output Parameter:
1017: . numDof - the number of dof
1019: Level: intermediate
1021: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
1022: @*/
1023: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1024: {
1025: PetscFunctionBegin;
1027: PetscAssertPointer(numDof, 4);
1028: PetscSectionCheckValidField(field, s->numFields);
1029: PetscCall(PetscSectionGetDof(s->field[field], point, numDof));
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: /*@
1034: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
1036: Not Collective
1038: Input Parameters:
1039: + s - the `PetscSection`
1040: . point - the point
1041: . field - the field
1042: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
1044: Level: intermediate
1046: Note:
1047: 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
1048: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1050: This is equivalent to
1051: .vb
1052: PetscSection fs;
1053: PetscSectionGetField(s,field,&fs)
1054: PetscSectionSetDof(fs,numDof)
1055: .ve
1057: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`, `PetscSectionAddDof()`, `PetscSectionSetDof()`
1058: @*/
1059: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1060: {
1061: PetscFunctionBegin;
1063: PetscSectionCheckValidField(field, s->numFields);
1064: PetscCall(PetscSectionSetDof(s->field[field], point, numDof));
1065: PetscFunctionReturn(PETSC_SUCCESS);
1066: }
1068: /*@
1069: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
1071: Not Collective
1073: Input Parameters:
1074: + s - the `PetscSection`
1075: . point - the point
1076: . field - the field
1077: - numDof - the number of dof
1079: Level: intermediate
1081: Notes:
1082: 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
1083: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1085: This is equivalent to
1086: .vb
1087: PetscSection fs;
1088: PetscSectionGetField(s,field,&fs)
1089: PetscSectionAddDof(fs,numDof)
1090: .ve
1092: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
1093: @*/
1094: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1095: {
1096: PetscFunctionBegin;
1098: PetscSectionCheckValidField(field, s->numFields);
1099: PetscCall(PetscSectionAddDof(s->field[field], point, numDof));
1100: PetscFunctionReturn(PETSC_SUCCESS);
1101: }
1103: /*@
1104: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
1106: Not Collective
1108: Input Parameters:
1109: + s - the `PetscSection`
1110: - point - the point
1112: Output Parameter:
1113: . numDof - the number of dof which are fixed by constraints
1115: Level: intermediate
1117: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
1118: @*/
1119: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
1120: {
1121: PetscFunctionBegin;
1123: PetscAssertPointer(numDof, 3);
1124: if (s->bc) PetscCall(PetscSectionGetDof(s->bc, point, numDof));
1125: else *numDof = 0;
1126: PetscFunctionReturn(PETSC_SUCCESS);
1127: }
1129: /*@
1130: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
1132: Not Collective
1134: Input Parameters:
1135: + s - the `PetscSection`
1136: . point - the point
1137: - numDof - the number of dof which are fixed by constraints
1139: Level: intermediate
1141: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1142: @*/
1143: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1144: {
1145: PetscFunctionBegin;
1147: if (numDof) {
1148: PetscCall(PetscSectionCheckConstraints_Private(s));
1149: PetscCall(PetscSectionSetDof(s->bc, point, numDof));
1150: }
1151: PetscFunctionReturn(PETSC_SUCCESS);
1152: }
1154: /*@
1155: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
1157: Not Collective
1159: Input Parameters:
1160: + s - the `PetscSection`
1161: . point - the point
1162: - numDof - the number of additional dof which are fixed by constraints
1164: Level: intermediate
1166: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1167: @*/
1168: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1169: {
1170: PetscFunctionBegin;
1172: if (numDof) {
1173: PetscCall(PetscSectionCheckConstraints_Private(s));
1174: PetscCall(PetscSectionAddDof(s->bc, point, numDof));
1175: }
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: /*@
1180: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
1182: Not Collective
1184: Input Parameters:
1185: + s - the `PetscSection`
1186: . point - the point
1187: - field - the field
1189: Output Parameter:
1190: . numDof - the number of dof which are fixed by constraints
1192: Level: intermediate
1194: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1195: @*/
1196: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1197: {
1198: PetscFunctionBegin;
1200: PetscAssertPointer(numDof, 4);
1201: PetscSectionCheckValidField(field, s->numFields);
1202: PetscCall(PetscSectionGetConstraintDof(s->field[field], point, numDof));
1203: PetscFunctionReturn(PETSC_SUCCESS);
1204: }
1206: /*@
1207: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1209: Not Collective
1211: Input Parameters:
1212: + s - the `PetscSection`
1213: . point - the point
1214: . field - the field
1215: - numDof - the number of dof which are fixed by constraints
1217: Level: intermediate
1219: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1220: @*/
1221: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1222: {
1223: PetscFunctionBegin;
1225: PetscSectionCheckValidField(field, s->numFields);
1226: PetscCall(PetscSectionSetConstraintDof(s->field[field], point, numDof));
1227: PetscFunctionReturn(PETSC_SUCCESS);
1228: }
1230: /*@
1231: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1233: Not Collective
1235: Input Parameters:
1236: + s - the `PetscSection`
1237: . point - the point
1238: . field - the field
1239: - numDof - the number of additional dof which are fixed by constraints
1241: Level: intermediate
1243: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1244: @*/
1245: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1246: {
1247: PetscFunctionBegin;
1249: PetscSectionCheckValidField(field, s->numFields);
1250: PetscCall(PetscSectionAddConstraintDof(s->field[field], point, numDof));
1251: PetscFunctionReturn(PETSC_SUCCESS);
1252: }
1254: /*@
1255: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1257: Not Collective
1259: Input Parameter:
1260: . s - the `PetscSection`
1262: Level: advanced
1264: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1265: @*/
1266: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1267: {
1268: PetscFunctionBegin;
1270: if (s->bc) {
1271: const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;
1273: PetscCall(PetscSectionSetUp(s->bc));
1274: if (last >= 0) PetscCall(PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices));
1275: else s->bcIndices = NULL;
1276: }
1277: PetscFunctionReturn(PETSC_SUCCESS);
1278: }
1280: /*@
1281: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point in preparation for use of the `PetscSection`
1283: Not Collective
1285: Input Parameter:
1286: . s - the `PetscSection`
1288: Level: intermediate
1290: Notes:
1291: If used, `PetscSectionSetPermutation()` must be called before this routine.
1293: `PetscSectionSetPointMajor()`, cannot be called after this routine.
1295: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
1296: @*/
1297: PetscErrorCode PetscSectionSetUp(PetscSection s)
1298: {
1299: PetscInt f;
1300: const PetscInt *pind = NULL;
1301: PetscCount offset = 0;
1303: PetscFunctionBegin;
1305: if (s->setup) PetscFunctionReturn(PETSC_SUCCESS);
1306: s->setup = PETSC_TRUE;
1307: /* Set offsets and field offsets for all points */
1308: /* Assume that all fields have the same chart */
1309: PetscCheck(s->includesConstraints, PETSC_COMM_SELF, PETSC_ERR_SUP, "PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE");
1310: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1311: if (s->pointMajor) {
1312: PetscCount foff;
1313: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1314: const PetscInt q = pind ? pind[p] : p;
1316: /* Set point offset */
1317: PetscCall(PetscIntCast(offset, &s->atlasOff[q]));
1318: offset += s->atlasDof[q];
1319: /* Set field offset */
1320: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1321: PetscSection sf = s->field[f];
1323: PetscCall(PetscIntCast(foff, &sf->atlasOff[q]));
1324: foff += sf->atlasDof[q];
1325: }
1326: }
1327: } else {
1328: /* Set field offsets for all points */
1329: for (f = 0; f < s->numFields; ++f) {
1330: PetscSection sf = s->field[f];
1332: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1333: const PetscInt q = pind ? pind[p] : p;
1335: PetscCall(PetscIntCast(offset, &sf->atlasOff[q]));
1336: offset += sf->atlasDof[q];
1337: }
1338: }
1339: /* Disable point offsets since these are unused */
1340: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1341: }
1342: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1343: /* Setup BC sections */
1344: PetscCall(PetscSectionSetUpBC(s));
1345: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetUpBC(s->field[f]));
1346: PetscFunctionReturn(PETSC_SUCCESS);
1347: }
1349: /*@
1350: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the `PetscSection`
1352: Not Collective
1354: Input Parameter:
1355: . s - the `PetscSection`
1357: Output Parameter:
1358: . maxDof - the maximum dof
1360: Level: intermediate
1362: Notes:
1363: The returned number is up-to-date without need for `PetscSectionSetUp()`.
1365: 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
1366: the maximum over all points of the value returned by `PetscSectionGetDof()` on this MPI process
1368: Developer Notes:
1369: The returned number is calculated lazily and stashed.
1371: A call to `PetscSectionInvalidateMaxDof_Internal()` invalidates the stashed value.
1373: `PetscSectionInvalidateMaxDof_Internal()` is called in `PetscSectionSetDof()`, `PetscSectionAddDof()` and `PetscSectionReset()`
1375: It should also be called every time `atlasDof` is modified directly.
1377: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1378: @*/
1379: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1380: {
1381: PetscInt p;
1383: PetscFunctionBegin;
1385: PetscAssertPointer(maxDof, 2);
1386: if (s->maxDof == PETSC_INT_MIN) {
1387: s->maxDof = 0;
1388: for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1389: }
1390: *maxDof = s->maxDof;
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1394: /*@
1395: PetscSectionGetStorageSize - Return the size of an array or local `Vec` capable of holding all the degrees of freedom defined in a `PetscSection`
1397: Not Collective
1399: Input Parameter:
1400: . s - the `PetscSection`
1402: Output Parameter:
1403: . size - the size of an array which can hold all the dofs
1405: Level: intermediate
1407: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1408: @*/
1409: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1410: {
1411: PetscInt64 n = 0;
1413: PetscFunctionBegin;
1415: PetscAssertPointer(size, 2);
1416: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1417: PetscCall(PetscIntCast(n, size));
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: }
1421: /*@
1422: PetscSectionGetConstrainedStorageSize - Return the size of an array or local `Vec` capable of holding all unconstrained degrees of freedom in a `PetscSection`
1424: Not Collective
1426: Input Parameter:
1427: . s - the `PetscSection`
1429: Output Parameter:
1430: . size - the size of an array which can hold all unconstrained dofs
1432: Level: intermediate
1434: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1435: @*/
1436: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1437: {
1438: PetscInt64 n = 0;
1440: PetscFunctionBegin;
1442: PetscAssertPointer(size, 2);
1443: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1444: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1445: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1446: }
1447: PetscCall(PetscIntCast(n, size));
1448: PetscFunctionReturn(PETSC_SUCCESS);
1449: }
1451: /*@
1452: PetscSectionCreateGlobalSection - Create a parallel section describing the global layout using
1453: a local (sequential) `PetscSection` on each MPI process and a `PetscSF` describing the section point overlap.
1455: Input Parameters:
1456: + s - The `PetscSection` for the local field layout
1457: . sf - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1458: . usePermutation - By default this is `PETSC_TRUE`, meaning any permutation of the local section is transferred to the global section
1459: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1460: - localOffsets - If `PETSC_TRUE`, use local rather than global offsets for the points
1462: Output Parameter:
1463: . gsection - The `PetscSection` for the global field layout
1465: Level: intermediate
1467: Notes:
1468: On each MPI process `gsection` inherits the chart of the `s` on that process.
1470: 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`.
1471: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1473: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCreateGlobalSectionCensored()`
1474: @*/
1475: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool usePermutation, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1476: {
1477: PetscSection gs;
1478: const PetscInt *pind = NULL;
1479: PetscInt *recv = NULL, *neg = NULL;
1480: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1481: PetscInt numFields, f, numComponents;
1482: PetscInt foff;
1484: PetscFunctionBegin;
1490: PetscAssertPointer(gsection, 6);
1491: PetscCheck(s->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
1492: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs));
1493: PetscCall(PetscSectionGetNumFields(s, &numFields));
1494: if (numFields > 0) PetscCall(PetscSectionSetNumFields(gs, numFields));
1495: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1496: PetscCall(PetscSectionSetChart(gs, pStart, pEnd));
1497: gs->includesConstraints = includeConstraints;
1498: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1499: nlocal = nroots; /* The local/leaf space matches global/root space */
1500: /* Must allocate for all points visible to SF, which may be more than this section */
1501: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1502: PetscCall(PetscSFGetLeafRange(sf, NULL, &maxleaf));
1503: PetscCheck(nroots >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %" PetscInt_FMT " < pEnd %" PetscInt_FMT, nroots, pEnd);
1504: PetscCheck(maxleaf < nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %" PetscInt_FMT " >= nroots %" PetscInt_FMT, maxleaf, nroots);
1505: PetscCall(PetscMalloc2(nroots, &neg, nlocal, &recv));
1506: PetscCall(PetscArrayzero(neg, nroots));
1507: }
1508: /* Mark all local points with negative dof */
1509: for (p = pStart; p < pEnd; ++p) {
1510: PetscCall(PetscSectionGetDof(s, p, &dof));
1511: PetscCall(PetscSectionSetDof(gs, p, dof));
1512: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1513: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(gs, p, cdof));
1514: if (neg) neg[p] = -(dof + 1);
1515: }
1516: PetscCall(PetscSectionSetUpBC(gs));
1517: 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]));
1518: if (nroots >= 0) {
1519: PetscCall(PetscArrayzero(recv, nlocal));
1520: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1521: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1522: for (p = pStart; p < pEnd; ++p) {
1523: if (recv[p] < 0) {
1524: gs->atlasDof[p - pStart] = recv[p];
1525: PetscCall(PetscSectionGetDof(s, p, &dof));
1526: 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);
1527: }
1528: }
1529: }
1530: /* Calculate new sizes, get process offset, and calculate point offsets */
1531: if (usePermutation && s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1532: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1533: const PetscInt q = pind ? pind[p] : p;
1535: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1536: gs->atlasOff[q] = off;
1537: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1538: }
1539: if (!localOffsets) {
1540: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)sf)));
1541: globalOff -= off;
1542: }
1543: for (p = pStart, off = 0; p < pEnd; ++p) {
1544: gs->atlasOff[p - pStart] += globalOff;
1545: if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1546: }
1547: if (usePermutation && s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1548: /* Put in negative offsets for ghost points */
1549: if (nroots >= 0) {
1550: PetscCall(PetscArrayzero(recv, nlocal));
1551: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1552: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1553: for (p = pStart; p < pEnd; ++p) {
1554: if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1555: }
1556: }
1557: PetscCall(PetscFree2(neg, recv));
1558: /* Set field dofs/offsets/constraints */
1559: for (f = 0; f < numFields; ++f) {
1560: const char *name;
1562: gs->field[f]->includesConstraints = includeConstraints;
1563: PetscCall(PetscSectionGetFieldComponents(s, f, &numComponents));
1564: PetscCall(PetscSectionSetFieldComponents(gs, f, numComponents));
1565: PetscCall(PetscSectionGetFieldName(s, f, &name));
1566: PetscCall(PetscSectionSetFieldName(gs, f, name));
1567: }
1568: for (p = pStart; p < pEnd; ++p) {
1569: PetscCall(PetscSectionGetOffset(gs, p, &off));
1570: for (f = 0, foff = off; f < numFields; ++f) {
1571: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
1572: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetFieldConstraintDof(gs, p, f, cdof));
1573: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
1574: PetscCall(PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof));
1575: PetscCall(PetscSectionSetFieldOffset(gs, p, f, foff));
1576: PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof));
1577: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1578: }
1579: }
1580: for (f = 0; f < numFields; ++f) {
1581: PetscSection gfs = gs->field[f];
1583: PetscCall(PetscSectionSetUpBC(gfs));
1584: 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]));
1585: }
1586: gs->setup = PETSC_TRUE;
1587: PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view"));
1588: *gsection = gs;
1589: PetscFunctionReturn(PETSC_SUCCESS);
1590: }
1592: /*@
1593: PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the globallayout using
1594: a local (sequential) `PetscSection` on each MPI process and an `PetscSF` describing the section point overlap.
1596: Input Parameters:
1597: + s - The `PetscSection` for the local field layout
1598: . sf - The `PetscSF` describing parallel layout of the section points
1599: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global vector will not possess constrained dofs
1600: . numExcludes - The number of exclusion ranges, this must have the same value on all MPI processes
1601: - 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
1603: Output Parameter:
1604: . gsection - The `PetscSection` for the global field layout
1606: Level: advanced
1608: Notes:
1609: On each MPI process `gsection` inherits the chart of the `s` on that process.
1611: 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`.
1612: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1614: This routine augments `PetscSectionCreateGlobalSection()` by allowing one to exclude certain ranges in the chart of the `PetscSection`
1616: Developer Notes:
1617: This is a terrible function name
1619: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`
1620: @*/
1621: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1622: {
1623: const PetscInt *pind = NULL;
1624: PetscInt *neg = NULL, *tmpOff = NULL;
1625: PetscInt pStart, pEnd, p, e, dof, cdof, globalOff = 0, nroots;
1626: PetscInt off;
1628: PetscFunctionBegin;
1631: PetscAssertPointer(gsection, 6);
1632: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
1633: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1634: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
1635: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1636: if (nroots >= 0) {
1637: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
1638: PetscCall(PetscCalloc1(nroots, &neg));
1639: if (nroots > pEnd - pStart) {
1640: PetscCall(PetscCalloc1(nroots, &tmpOff));
1641: } else {
1642: tmpOff = &(*gsection)->atlasDof[-pStart];
1643: }
1644: }
1645: /* Mark ghost points with negative dof */
1646: for (p = pStart; p < pEnd; ++p) {
1647: for (e = 0; e < numExcludes; ++e) {
1648: if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1649: PetscCall(PetscSectionSetDof(*gsection, p, 0));
1650: break;
1651: }
1652: }
1653: if (e < numExcludes) continue;
1654: PetscCall(PetscSectionGetDof(s, p, &dof));
1655: PetscCall(PetscSectionSetDof(*gsection, p, dof));
1656: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1657: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
1658: if (neg) neg[p] = -(dof + 1);
1659: }
1660: PetscCall(PetscSectionSetUpBC(*gsection));
1661: if (nroots >= 0) {
1662: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1663: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1664: if (nroots > pEnd - pStart) {
1665: for (p = pStart; p < pEnd; ++p) {
1666: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1667: }
1668: }
1669: }
1670: /* Calculate new sizes, get process offset, and calculate point offsets */
1671: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1672: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1673: const PetscInt q = pind ? pind[p] : p;
1675: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1676: (*gsection)->atlasOff[q] = off;
1677: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1678: }
1679: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1680: globalOff -= off;
1681: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1682: (*gsection)->atlasOff[p] += globalOff;
1683: if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1684: }
1685: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1686: /* Put in negative offsets for ghost points */
1687: if (nroots >= 0) {
1688: if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1689: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1690: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1691: if (nroots > pEnd - pStart) {
1692: for (p = pStart; p < pEnd; ++p) {
1693: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1694: }
1695: }
1696: }
1697: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
1698: PetscCall(PetscFree(neg));
1699: PetscFunctionReturn(PETSC_SUCCESS);
1700: }
1702: /*@
1703: PetscSectionGetPointLayout - Get a `PetscLayout` for the points with nonzero dof counts of the unnamed default field within this `PetscSection`s local chart
1705: Collective
1707: Input Parameters:
1708: + comm - The `MPI_Comm`
1709: - s - The `PetscSection`
1711: Output Parameter:
1712: . layout - The point layout for the data that defines the section
1714: Level: advanced
1716: Notes:
1717: `PetscSectionGetValueLayout()` provides similar information but counting the total number of degrees of freedom on the MPI process (excluding constrained
1718: degrees of freedom).
1720: This count includes constrained degrees of freedom
1722: This is usually called on the default global section.
1724: Example:
1725: .vb
1726: The chart is [2,5), point 2 has 2 dof, point 3 has 0 dof, point 4 has 1 dof
1727: The local size of the `PetscLayout` is 2 since 2 points have a non-zero number of dof
1728: .ve
1730: Developer Notes:
1731: I find the names of these two functions extremely non-informative
1733: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1734: @*/
1735: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1736: {
1737: PetscInt pStart, pEnd, p, localSize = 0;
1739: PetscFunctionBegin;
1740: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1741: for (p = pStart; p < pEnd; ++p) {
1742: PetscInt dof;
1744: PetscCall(PetscSectionGetDof(s, p, &dof));
1745: if (dof >= 0) ++localSize;
1746: }
1747: PetscCall(PetscLayoutCreate(comm, layout));
1748: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1749: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1750: PetscCall(PetscLayoutSetUp(*layout));
1751: PetscFunctionReturn(PETSC_SUCCESS);
1752: }
1754: /*@
1755: PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs of a `PetscSection`
1757: Collective
1759: Input Parameters:
1760: + comm - The `MPI_Comm`
1761: - s - The `PetscSection`
1763: Output Parameter:
1764: . layout - The dof layout for the section
1766: Level: advanced
1768: Notes:
1769: `PetscSectionGetPointLayout()` provides similar information but only counting the number of points with nonzero degrees of freedom and
1770: including the constrained degrees of freedom
1772: This is usually called for the default global section.
1774: Example:
1775: .vb
1776: The chart is [2,5), point 2 has 4 dof (2 constrained), point 3 has 0 dof, point 4 has 1 dof (not constrained)
1777: The local size of the `PetscLayout` is 3 since there are 3 unconstrained degrees of freedom on this MPI process
1778: .ve
1780: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1781: @*/
1782: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1783: {
1784: PetscInt pStart, pEnd, p, localSize = 0;
1786: PetscFunctionBegin;
1788: PetscAssertPointer(layout, 3);
1789: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1790: for (p = pStart; p < pEnd; ++p) {
1791: PetscInt dof, cdof;
1793: PetscCall(PetscSectionGetDof(s, p, &dof));
1794: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1795: if (dof - cdof > 0) localSize += dof - cdof;
1796: }
1797: PetscCall(PetscLayoutCreate(comm, layout));
1798: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1799: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1800: PetscCall(PetscLayoutSetUp(*layout));
1801: PetscFunctionReturn(PETSC_SUCCESS);
1802: }
1804: /*@
1805: PetscSectionGetOffset - Return the offset into an array or `Vec` for the dof associated with the given point.
1807: Not Collective
1809: Input Parameters:
1810: + s - the `PetscSection`
1811: - point - the point
1813: Output Parameter:
1814: . offset - the offset
1816: Level: intermediate
1818: Notes:
1819: In a global section, `offset` will be negative for points not owned by this process.
1821: This is for the unnamed default field in the `PetscSection` not the named fields
1823: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1825: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetPointMajor()`
1826: @*/
1827: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1828: {
1829: PetscFunctionBegin;
1831: PetscAssertPointer(offset, 3);
1832: 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);
1833: *offset = s->atlasOff[point - s->pStart];
1834: PetscFunctionReturn(PETSC_SUCCESS);
1835: }
1837: /*@
1838: PetscSectionSetOffset - Set the offset into an array or `Vec` for the dof associated with the given point.
1840: Not Collective
1842: Input Parameters:
1843: + s - the `PetscSection`
1844: . point - the point
1845: - offset - the offset, these values may be negative indicating the values are off process
1847: Level: developer
1849: Note:
1850: The user usually does not call this function, but uses `PetscSectionSetUp()`
1852: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1853: @*/
1854: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1855: {
1856: PetscFunctionBegin;
1858: 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);
1859: s->atlasOff[point - s->pStart] = offset;
1860: PetscFunctionReturn(PETSC_SUCCESS);
1861: }
1863: /*@
1864: PetscSectionGetFieldOffset - Return the offset into an array or `Vec` for the field dof associated with the given point.
1866: Not Collective
1868: Input Parameters:
1869: + s - the `PetscSection`
1870: . point - the point
1871: - field - the field
1873: Output Parameter:
1874: . offset - the offset
1876: Level: intermediate
1878: Notes:
1879: In a global section, `offset` will be negative for points not owned by this process.
1881: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1883: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldPointOffset()`
1884: @*/
1885: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1886: {
1887: PetscFunctionBegin;
1889: PetscAssertPointer(offset, 4);
1890: PetscSectionCheckValidField(field, s->numFields);
1891: PetscCall(PetscSectionGetOffset(s->field[field], point, offset));
1892: PetscFunctionReturn(PETSC_SUCCESS);
1893: }
1895: /*@
1896: PetscSectionSetFieldOffset - Set the offset into an array or `Vec` for the dof associated with the given field at a point.
1898: Not Collective
1900: Input Parameters:
1901: + s - the `PetscSection`
1902: . point - the point
1903: . field - the field
1904: - offset - the offset, these values may be negative indicating the values are off process
1906: Level: developer
1908: Note:
1909: The user usually does not call this function, but uses `PetscSectionSetUp()`
1911: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1912: @*/
1913: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1914: {
1915: PetscFunctionBegin;
1917: PetscSectionCheckValidField(field, s->numFields);
1918: PetscCall(PetscSectionSetOffset(s->field[field], point, offset));
1919: PetscFunctionReturn(PETSC_SUCCESS);
1920: }
1922: /*@
1923: PetscSectionGetFieldPointOffset - Return the offset for the first field dof associated with the given point relative to the offset for that point for the
1924: unnamed default field's first dof
1926: Not Collective
1928: Input Parameters:
1929: + s - the `PetscSection`
1930: . point - the point
1931: - field - the field
1933: Output Parameter:
1934: . offset - the offset
1936: Level: advanced
1938: Note:
1939: This ignores constraints
1941: Example:
1942: .vb
1943: if PetscSectionSetPointMajor(s,PETSC_TRUE)
1944: The unnamed default field has 3 dof at `point`
1945: Field 0 has 2 dof at `point`
1946: Then PetscSectionGetFieldPointOffset(s,point,1,&offset) returns and offset of 5
1947: .ve
1949: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldOffset()`
1950: @*/
1951: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1952: {
1953: PetscInt off, foff;
1955: PetscFunctionBegin;
1957: PetscAssertPointer(offset, 4);
1958: PetscSectionCheckValidField(field, s->numFields);
1959: PetscCall(PetscSectionGetOffset(s, point, &off));
1960: PetscCall(PetscSectionGetOffset(s->field[field], point, &foff));
1961: *offset = foff - off;
1962: PetscFunctionReturn(PETSC_SUCCESS);
1963: }
1965: /*@
1966: PetscSectionGetOffsetRange - Return the full range of offsets [`start`, `end`) for a `PetscSection`
1968: Not Collective
1970: Input Parameter:
1971: . s - the `PetscSection`
1973: Output Parameters:
1974: + start - the minimum offset
1975: - end - one more than the maximum offset
1977: Level: intermediate
1979: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1980: @*/
1981: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1982: {
1983: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1985: PetscFunctionBegin;
1987: if (s->atlasOff) {
1988: os = s->atlasOff[0];
1989: oe = s->atlasOff[0];
1990: }
1991: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1992: for (p = 0; p < pEnd - pStart; ++p) {
1993: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1995: if (off >= 0) {
1996: os = PetscMin(os, off);
1997: oe = PetscMax(oe, off + dof);
1998: }
1999: }
2000: if (start) *start = os;
2001: if (end) *end = oe;
2002: PetscFunctionReturn(PETSC_SUCCESS);
2003: }
2005: /*@
2006: PetscSectionCreateSubsection - Create a new, smaller `PetscSection` composed of only selected fields
2008: Collective
2010: Input Parameters:
2011: + s - the `PetscSection`
2012: . len - the number of subfields
2013: - fields - the subfield numbers
2015: Output Parameter:
2016: . subs - the subsection
2018: Level: advanced
2020: Notes:
2021: The chart of `subs` is the same as the chart of `s`
2023: This will error if a fieldnumber is out of range
2025: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2026: @*/
2027: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
2028: {
2029: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
2031: PetscFunctionBegin;
2032: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2034: PetscAssertPointer(fields, 3);
2035: PetscAssertPointer(subs, 4);
2036: PetscCall(PetscSectionGetNumFields(s, &nF));
2037: 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);
2038: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2039: PetscCall(PetscSectionSetNumFields(*subs, len));
2040: for (f = 0; f < len; ++f) {
2041: const char *name = NULL;
2042: PetscInt numComp = 0;
2043: PetscSectionSym sym;
2045: PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
2046: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2047: PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
2048: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2049: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
2050: PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
2051: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2052: }
2053: PetscCall(PetscSectionGetFieldSym(s, fields[f], &sym));
2054: PetscCall(PetscSectionSetFieldSym(*subs, f, sym));
2055: }
2056: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2057: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2058: for (p = pStart; p < pEnd; ++p) {
2059: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
2061: for (f = 0; f < len; ++f) {
2062: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2063: PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
2064: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2065: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
2066: dof += fdof;
2067: cdof += cfdof;
2068: }
2069: PetscCall(PetscSectionSetDof(*subs, p, dof));
2070: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
2071: maxCdof = PetscMax(cdof, maxCdof);
2072: }
2073: PetscBT bst, subbst;
2075: PetscCall(PetscSectionGetBlockStarts(s, &bst));
2076: if (bst) {
2077: PetscCall(PetscBTCreate(pEnd - pStart, &subbst));
2078: PetscCall(PetscBTCopy(subbst, pEnd - pStart, bst));
2079: PetscCall(PetscSectionSetBlockStarts(*subs, subbst));
2080: }
2081: PetscCall(PetscSectionSetUp(*subs));
2082: if (maxCdof) {
2083: PetscInt *indices;
2085: PetscCall(PetscMalloc1(maxCdof, &indices));
2086: for (p = pStart; p < pEnd; ++p) {
2087: PetscInt cdof;
2089: PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
2090: if (cdof) {
2091: const PetscInt *oldIndices = NULL;
2092: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
2094: for (f = 0; f < len; ++f) {
2095: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2096: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2097: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
2098: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
2099: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
2100: numConst += cfdof;
2101: fOff += fdof;
2102: }
2103: PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
2104: }
2105: }
2106: PetscCall(PetscFree(indices));
2107: }
2108: PetscFunctionReturn(PETSC_SUCCESS);
2109: }
2111: /*@
2112: PetscSectionCreateComponentSubsection - Create a new, smaller `PetscSection` composed of only selected components
2114: Collective
2116: Input Parameters:
2117: + s - the `PetscSection`
2118: . len - the number of components
2119: - comps - the component numbers
2121: Output Parameter:
2122: . subs - the subsection
2124: Level: advanced
2126: Notes:
2127: The chart of `subs` is the same as the chart of `s`
2129: This will error if the section has more than one field, or if a component number is out of range
2131: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2132: @*/
2133: PetscErrorCode PetscSectionCreateComponentSubsection(PetscSection s, PetscInt len, const PetscInt comps[], PetscSection *subs)
2134: {
2135: PetscSectionSym sym;
2136: const char *name = NULL;
2137: PetscInt Nf, pStart, pEnd;
2139: PetscFunctionBegin;
2140: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2142: PetscAssertPointer(comps, 3);
2143: PetscAssertPointer(subs, 4);
2144: PetscCall(PetscSectionGetNumFields(s, &Nf));
2145: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "This method can only handle one field, not %" PetscInt_FMT, Nf);
2146: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2147: PetscCall(PetscSectionSetNumFields(*subs, 1));
2148: PetscCall(PetscSectionGetFieldName(s, 0, &name));
2149: PetscCall(PetscSectionSetFieldName(*subs, 0, name));
2150: PetscCall(PetscSectionSetFieldComponents(*subs, 0, len));
2151: PetscCall(PetscSectionGetFieldSym(s, 0, &sym));
2152: PetscCall(PetscSectionSetFieldSym(*subs, 0, sym));
2153: for (PetscInt c = 0; c < len; ++c) {
2154: PetscCall(PetscSectionGetComponentName(s, 0, comps[c], &name));
2155: PetscCall(PetscSectionSetComponentName(*subs, 0, c, name));
2156: }
2157: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2158: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2159: for (PetscInt p = pStart; p < pEnd; ++p) {
2160: PetscInt dof, cdof, cfdof;
2162: PetscCall(PetscSectionGetDof(s, p, &dof));
2163: if (!dof) continue;
2164: PetscCall(PetscSectionGetFieldConstraintDof(s, p, 0, &cfdof));
2165: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2166: PetscCheck(!cdof && !cfdof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Component selection does not work with constraints");
2167: PetscCall(PetscSectionSetFieldDof(*subs, p, 0, len));
2168: PetscCall(PetscSectionSetDof(*subs, p, len));
2169: }
2170: PetscCall(PetscSectionSetUp(*subs));
2171: PetscFunctionReturn(PETSC_SUCCESS);
2172: }
2174: /*@
2175: PetscSectionCreateSupersection - Create a new, larger section composed of multiple `PetscSection`s
2177: Collective
2179: Input Parameters:
2180: + s - the input sections
2181: - len - the number of input sections
2183: Output Parameter:
2184: . supers - the supersection
2186: Level: advanced
2188: Notes:
2189: The section offsets now refer to a new, larger vector.
2191: Developer Notes:
2192: Needs to explain how the sections are composed
2194: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2195: @*/
2196: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2197: {
2198: PetscInt Nf = 0, f, pStart = PETSC_INT_MAX, pEnd = 0, p, maxCdof = 0, i;
2200: PetscFunctionBegin;
2201: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2202: for (i = 0; i < len; ++i) {
2203: PetscInt nf, pStarti, pEndi;
2205: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2206: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2207: pStart = PetscMin(pStart, pStarti);
2208: pEnd = PetscMax(pEnd, pEndi);
2209: Nf += nf;
2210: }
2211: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2212: PetscCall(PetscSectionSetNumFields(*supers, Nf));
2213: for (i = 0, f = 0; i < len; ++i) {
2214: PetscInt nf, fi, ci;
2216: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2217: for (fi = 0; fi < nf; ++fi, ++f) {
2218: const char *name = NULL;
2219: PetscInt numComp = 0;
2221: PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2222: PetscCall(PetscSectionSetFieldName(*supers, f, name));
2223: PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2224: PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2225: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2226: PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2227: PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2228: }
2229: }
2230: }
2231: PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2232: for (p = pStart; p < pEnd; ++p) {
2233: PetscInt dof = 0, cdof = 0;
2235: for (i = 0, f = 0; i < len; ++i) {
2236: PetscInt nf, fi, pStarti, pEndi;
2237: PetscInt fdof = 0, cfdof = 0;
2239: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2240: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2241: if ((p < pStarti) || (p >= pEndi)) continue;
2242: for (fi = 0; fi < nf; ++fi, ++f) {
2243: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2244: PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2245: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2246: if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2247: dof += fdof;
2248: cdof += cfdof;
2249: }
2250: }
2251: PetscCall(PetscSectionSetDof(*supers, p, dof));
2252: if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2253: maxCdof = PetscMax(cdof, maxCdof);
2254: }
2255: PetscCall(PetscSectionSetUp(*supers));
2256: if (maxCdof) {
2257: PetscInt *indices;
2259: PetscCall(PetscMalloc1(maxCdof, &indices));
2260: for (p = pStart; p < pEnd; ++p) {
2261: PetscInt cdof;
2263: PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2264: if (cdof) {
2265: PetscInt dof, numConst = 0, fOff = 0;
2267: for (i = 0, f = 0; i < len; ++i) {
2268: const PetscInt *oldIndices = NULL;
2269: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
2271: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2272: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2273: if ((p < pStarti) || (p >= pEndi)) continue;
2274: for (fi = 0; fi < nf; ++fi, ++f) {
2275: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2276: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2277: PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2278: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2279: PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2280: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2281: numConst += cfdof;
2282: }
2283: PetscCall(PetscSectionGetDof(s[i], p, &dof));
2284: fOff += dof;
2285: }
2286: PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2287: }
2288: }
2289: PetscCall(PetscFree(indices));
2290: }
2291: PetscFunctionReturn(PETSC_SUCCESS);
2292: }
2294: static PetscErrorCode PetscSectionCreateSubplexSection_Private(PetscSection s, IS subpointIS, PetscBool renumberPoints, PetscSection *subs)
2295: {
2296: const PetscInt *points = NULL, *indices = NULL;
2297: PetscInt *spoints = NULL, *order = NULL;
2298: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
2300: PetscFunctionBegin;
2303: PetscAssertPointer(subs, 4);
2304: PetscCall(PetscSectionGetNumFields(s, &numFields));
2305: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2306: if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2307: for (f = 0; f < numFields; ++f) {
2308: const char *name = NULL;
2309: PetscInt numComp = 0;
2311: PetscCall(PetscSectionGetFieldName(s, f, &name));
2312: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2313: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2314: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2315: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2316: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2317: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2318: }
2319: }
2320: /* For right now, we do not try to squeeze the subchart */
2321: if (subpointIS) {
2322: PetscCall(ISGetLocalSize(subpointIS, &numSubpoints));
2323: PetscCall(ISGetIndices(subpointIS, &points));
2324: }
2325: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2326: if (renumberPoints) {
2327: PetscBool sorted;
2329: spStart = 0;
2330: spEnd = numSubpoints;
2331: PetscCall(ISSorted(subpointIS, &sorted));
2332: if (!sorted) {
2333: PetscCall(PetscMalloc2(numSubpoints, &spoints, numSubpoints, &order));
2334: PetscCall(PetscArraycpy(spoints, points, numSubpoints));
2335: for (PetscInt i = 0; i < numSubpoints; ++i) order[i] = i;
2336: PetscCall(PetscSortIntWithArray(numSubpoints, spoints, order));
2337: }
2338: } else {
2339: PetscCall(ISGetMinMax(subpointIS, &spStart, &spEnd));
2340: ++spEnd;
2341: }
2342: PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2343: for (p = pStart; p < pEnd; ++p) {
2344: PetscInt dof, cdof, fdof = 0, cfdof = 0;
2346: PetscCall(PetscFindInt(p, numSubpoints, spoints ? spoints : points, &subp));
2347: if (subp < 0) continue;
2348: if (!renumberPoints) subp = p;
2349: else subp = order ? order[subp] : subp;
2350: for (f = 0; f < numFields; ++f) {
2351: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2352: PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2353: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2354: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2355: }
2356: PetscCall(PetscSectionGetDof(s, p, &dof));
2357: PetscCall(PetscSectionSetDof(*subs, subp, dof));
2358: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2359: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2360: }
2361: PetscCall(PetscSectionSetUp(*subs));
2362: /* Change offsets to original offsets */
2363: for (p = pStart; p < pEnd; ++p) {
2364: PetscInt off, foff = 0;
2366: PetscCall(PetscFindInt(p, numSubpoints, spoints ? spoints : points, &subp));
2367: if (subp < 0) continue;
2368: if (!renumberPoints) subp = p;
2369: else subp = order ? order[subp] : subp;
2370: for (f = 0; f < numFields; ++f) {
2371: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2372: PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2373: }
2374: PetscCall(PetscSectionGetOffset(s, p, &off));
2375: PetscCall(PetscSectionSetOffset(*subs, subp, off));
2376: }
2377: /* Copy constraint indices */
2378: for (subp = spStart; subp < spEnd; ++subp) {
2379: PetscInt cdof;
2381: PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2382: if (cdof) {
2383: for (f = 0; f < numFields; ++f) {
2384: PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2385: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2386: }
2387: PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2388: PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2389: }
2390: }
2391: if (subpointIS) PetscCall(ISRestoreIndices(subpointIS, &points));
2392: PetscCall(PetscFree2(spoints, order));
2393: PetscFunctionReturn(PETSC_SUCCESS);
2394: }
2396: /*@
2397: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2399: Collective
2401: Input Parameters:
2402: + s - the `PetscSection`
2403: - subpointIS - a sorted list of points in the original mesh which are in the submesh
2405: Output Parameter:
2406: . subs - the subsection
2408: Level: advanced
2410: Notes:
2411: 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))`
2413: Compare this with `PetscSectionCreateSubdomainSection()` that does not map the points numbers to start at zero but leaves them as before
2415: Developer Notes:
2416: 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`
2418: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2419: @*/
2420: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointIS, PetscSection *subs)
2421: {
2422: PetscFunctionBegin;
2423: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointIS, PETSC_TRUE, subs));
2424: PetscFunctionReturn(PETSC_SUCCESS);
2425: }
2427: /*@
2428: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2430: Collective
2432: Input Parameters:
2433: + s - the `PetscSection`
2434: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2436: Output Parameter:
2437: . subs - the subsection
2439: Level: advanced
2441: Notes:
2442: 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`
2443: is `[min(subpointMap),max(subpointMap)+1)`
2445: Compare this with `PetscSectionCreateSubmeshSection()` that maps the point numbers to start at zero
2447: Developer Notes:
2448: 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`
2450: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2451: @*/
2452: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2453: {
2454: PetscFunctionBegin;
2455: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_FALSE, subs));
2456: PetscFunctionReturn(PETSC_SUCCESS);
2457: }
2459: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2460: {
2461: PetscInt p;
2462: PetscMPIInt rank;
2464: PetscFunctionBegin;
2465: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2466: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2467: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2468: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2469: if (s->bc && s->bc->atlasDof[p] > 0) {
2470: PetscInt b;
2471: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2472: if (s->bcIndices) {
2473: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2474: }
2475: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2476: } else {
2477: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2478: }
2479: }
2480: PetscCall(PetscViewerFlush(viewer));
2481: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2482: if (s->sym) {
2483: PetscCall(PetscViewerASCIIPushTab(viewer));
2484: PetscCall(PetscSectionSymView(s->sym, viewer));
2485: PetscCall(PetscViewerASCIIPopTab(viewer));
2486: }
2487: PetscFunctionReturn(PETSC_SUCCESS);
2488: }
2490: /*@
2491: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2493: Collective
2495: Input Parameters:
2496: + A - the `PetscSection` object to view
2497: . obj - Optional object that provides the options prefix used for the options
2498: - name - command line option
2500: Level: intermediate
2502: Note:
2503: See `PetscObjectViewFromOptions()` for available values of `PetscViewer` and `PetscViewerFormat`
2505: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2506: @*/
2507: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2508: {
2509: PetscFunctionBegin;
2511: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2512: PetscFunctionReturn(PETSC_SUCCESS);
2513: }
2515: /*@
2516: PetscSectionView - Views a `PetscSection`
2518: Collective
2520: Input Parameters:
2521: + s - the `PetscSection` object to view
2522: - viewer - the viewer
2524: Level: beginner
2526: Note:
2527: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2528: distribution independent data, such as dofs, offsets, constraint dofs,
2529: and constraint indices. Points that have negative dofs, for instance,
2530: are not saved as they represent points owned by other processes.
2531: Point numbering and rank assignment is currently not stored.
2532: The saved section can be loaded with `PetscSectionLoad()`.
2534: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2535: @*/
2536: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2537: {
2538: PetscBool isascii, ishdf5;
2539: PetscInt f;
2541: PetscFunctionBegin;
2543: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2545: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2546: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2547: if (isascii) {
2548: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2549: if (s->numFields) {
2550: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2551: for (f = 0; f < s->numFields; ++f) {
2552: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " \"%s\" with %" PetscInt_FMT " components\n", f, s->fieldNames[f], s->numFieldComponents[f]));
2553: PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2554: }
2555: } else {
2556: PetscCall(PetscSectionView_ASCII(s, viewer));
2557: }
2558: } else if (ishdf5) {
2559: #if PetscDefined(HAVE_HDF5)
2560: PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2561: #else
2562: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2563: #endif
2564: }
2565: PetscFunctionReturn(PETSC_SUCCESS);
2566: }
2568: /*@
2569: PetscSectionLoad - Loads a `PetscSection`
2571: Collective
2573: Input Parameters:
2574: + s - the `PetscSection` object to load
2575: - viewer - the viewer
2577: Level: beginner
2579: Note:
2580: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2581: a section saved with `PetscSectionView()`. The number of processes
2582: used here (N) does not need to be the same as that used when saving.
2583: After calling this function, the chart of s on rank i will be set
2584: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2585: saved section points.
2587: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2588: @*/
2589: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2590: {
2591: PetscBool ishdf5;
2593: PetscFunctionBegin;
2596: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2597: PetscCheck(ishdf5, PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2598: #if PetscDefined(HAVE_HDF5)
2599: PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2600: PetscFunctionReturn(PETSC_SUCCESS);
2601: #else
2602: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2603: #endif
2604: }
2606: static inline PetscErrorCode PrintArrayElement(void *array, PetscDataType data_type, PetscCount index, PetscViewer viewer)
2607: {
2608: PetscFunctionBeginUser;
2609: switch (data_type) {
2610: case PETSC_INT: {
2611: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt_FMT, ((PetscInt *)array)[index]));
2612: break;
2613: }
2614: case PETSC_INT32: {
2615: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt32_FMT, ((PetscInt32 *)array)[index]));
2616: break;
2617: }
2618: case PETSC_INT64: {
2619: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt64_FMT, ((PetscInt64 *)array)[index]));
2620: break;
2621: }
2622: case PETSC_COUNT: {
2623: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscCount_FMT, ((PetscCount *)array)[index]));
2624: break;
2625: }
2626: // PETSC_SCALAR is set to the appropriate type
2627: case PETSC_DOUBLE: {
2628: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", ((double *)array)[index]));
2629: break;
2630: }
2631: case PETSC_FLOAT: {
2632: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((float *)array)[index]));
2633: break;
2634: }
2635: #if defined(PETSC_USE_REAL___FLOAT128)
2636: case PETSC___FLOAT128: {
2637: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((PetscReal *)array)[index]));
2638: break;
2639: }
2640: #endif
2641: #if defined(PETSC_USE_REAL___FP16)
2642: case PETSC___FP16: {
2643: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((PetscReal *)array)[index]));
2644: break;
2645: }
2646: #endif
2647: #if defined(PETSC_HAVE_COMPLEX)
2648: case PETSC_COMPLEX: {
2649: PetscComplex v = ((PetscComplex *)array)[index];
2650: if (PetscImaginaryPartComplex(v) > 0.0) {
2651: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPartComplex(v), (double)PetscImaginaryPartComplex(v)));
2652: } else if (PetscImaginaryPartComplex(v) < 0.0) {
2653: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPartComplex(v), (double)(-PetscImaginaryPartComplex(v))));
2654: } else {
2655: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPartComplex(v)));
2656: }
2657: break;
2658: }
2659: #endif
2660: default:
2661: SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "PetscDataType %d (%s) not supported", data_type, PetscDataTypes[data_type]);
2662: }
2663: PetscFunctionReturn(PETSC_SUCCESS);
2664: }
2666: PetscErrorCode PetscSectionArrayView_ASCII_Internal(PetscSection s, void *array, PetscDataType data_type, PetscViewer viewer)
2667: {
2668: PetscInt p, i;
2669: PetscMPIInt rank;
2671: PetscFunctionBegin;
2672: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2673: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2674: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2675: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2676: if (s->bc && (s->bc->atlasDof[p] > 0)) {
2677: PetscInt b;
2679: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2680: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) PetscCall(PrintArrayElement(array, data_type, i, viewer));
2681: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " constrained"));
2682: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2683: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2684: } else {
2685: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2686: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) PetscCall(PrintArrayElement(array, data_type, i, viewer));
2687: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2688: }
2689: }
2690: PetscCall(PetscViewerFlush(viewer));
2691: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2692: PetscFunctionReturn(PETSC_SUCCESS);
2693: }
2695: /*@
2696: PetscSectionArrayView - View an array, using the section to structure the values
2698: Collective
2700: Input Parameters:
2701: + s - the organizing `PetscSection`
2702: . array - the array of values
2703: . data_type - the `PetscDataType` of the array
2704: - viewer - the `PetscViewer`
2706: Level: developer
2708: .seealso: `PetscSection`, `PetscViewer`, `PetscSectionCreate()`, `VecSetValuesSection()`, `PetscSectionVecView()`
2709: @*/
2710: PetscErrorCode PetscSectionArrayView(PetscSection s, void *array, PetscDataType data_type, PetscViewer viewer)
2711: {
2712: PetscBool isascii;
2713: PetscInt f;
2715: PetscFunctionBegin;
2717: if (!array) {
2718: PetscInt size;
2719: PetscCall(PetscSectionGetStorageSize(s, &size));
2720: PetscCheck(size == 0, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_SIZ, "NULL array passed, but section's storage size is non-zero");
2721: } else PetscAssertPointer(array, 2);
2722: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2724: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2725: if (isascii) {
2726: if (s->numFields) {
2727: PetscCall(PetscViewerASCIIPrintf(viewer, "Array with %" PetscInt_FMT " fields\n", s->numFields));
2728: for (f = 0; f < s->numFields; ++f) {
2729: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
2730: PetscCall(PetscSectionArrayView_ASCII_Internal(s->field[f], array, data_type, viewer));
2731: }
2732: } else {
2733: PetscCall(PetscSectionArrayView_ASCII_Internal(s, array, data_type, viewer));
2734: }
2735: }
2736: PetscFunctionReturn(PETSC_SUCCESS);
2737: }
2739: /*@
2740: PetscSectionResetClosurePermutation - Remove any existing closure permutation
2742: Input Parameter:
2743: . section - The `PetscSection`
2745: Level: intermediate
2747: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionSetClosureIndex()`, `PetscSectionReset()`
2748: @*/
2749: PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2750: {
2751: PetscSectionClosurePermVal clVal;
2753: PetscFunctionBegin;
2754: if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2755: kh_foreach_value(section->clHash, clVal, {
2756: PetscCall(PetscFree(clVal.perm));
2757: PetscCall(PetscFree(clVal.invPerm));
2758: });
2759: kh_destroy(ClPerm, section->clHash);
2760: section->clHash = NULL;
2761: PetscFunctionReturn(PETSC_SUCCESS);
2762: }
2764: /*@
2765: PetscSectionReset - Frees all section data, the section is then as if `PetscSectionCreate()` had just been called.
2767: Not Collective
2769: Input Parameter:
2770: . s - the `PetscSection`
2772: Level: beginner
2774: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`
2775: @*/
2776: PetscErrorCode PetscSectionReset(PetscSection s)
2777: {
2778: PetscInt f, c;
2780: PetscFunctionBegin;
2782: for (f = 0; f < s->numFields; ++f) {
2783: PetscCall(PetscSectionDestroy(&s->field[f]));
2784: PetscCall(PetscFree(s->fieldNames[f]));
2785: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2786: PetscCall(PetscFree(s->compNames[f]));
2787: }
2788: PetscCall(PetscFree(s->numFieldComponents));
2789: PetscCall(PetscFree(s->fieldNames));
2790: PetscCall(PetscFree(s->compNames));
2791: PetscCall(PetscFree(s->field));
2792: PetscCall(PetscSectionDestroy(&s->bc));
2793: PetscCall(PetscFree(s->bcIndices));
2794: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2795: PetscCall(PetscSectionDestroy(&s->clSection));
2796: PetscCall(ISDestroy(&s->clPoints));
2797: PetscCall(ISDestroy(&s->perm));
2798: PetscCall(PetscBTDestroy(&s->blockStarts));
2799: PetscCall(PetscSectionResetClosurePermutation(s));
2800: PetscCall(PetscSectionSymDestroy(&s->sym));
2801: PetscCall(PetscSectionDestroy(&s->clSection));
2802: PetscCall(ISDestroy(&s->clPoints));
2803: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2804: s->pStart = -1;
2805: s->pEnd = -1;
2806: s->maxDof = 0;
2807: s->setup = PETSC_FALSE;
2808: s->numFields = 0;
2809: s->clObj = NULL;
2810: PetscFunctionReturn(PETSC_SUCCESS);
2811: }
2813: /*@
2814: PetscSectionDestroy - Frees a `PetscSection`
2816: Not Collective
2818: Input Parameter:
2819: . s - the `PetscSection`
2821: Level: beginner
2823: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2824: @*/
2825: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2826: {
2827: PetscFunctionBegin;
2828: if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2830: if (--((PetscObject)*s)->refct > 0) {
2831: *s = NULL;
2832: PetscFunctionReturn(PETSC_SUCCESS);
2833: }
2834: PetscCall(PetscSectionReset(*s));
2835: PetscCall(PetscHeaderDestroy(s));
2836: PetscFunctionReturn(PETSC_SUCCESS);
2837: }
2839: static PetscErrorCode VecIntGetValuesSection_Private(const PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2840: {
2841: const PetscInt p = point - s->pStart;
2843: PetscFunctionBegin;
2845: *values = &baseArray[s->atlasOff[p]];
2846: PetscFunctionReturn(PETSC_SUCCESS);
2847: }
2849: static PetscErrorCode VecIntSetValuesSection_Private(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2850: {
2851: PetscInt *array;
2852: const PetscInt p = point - s->pStart;
2853: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2854: PetscInt cDim = 0;
2856: PetscFunctionBegin;
2858: PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2859: array = &baseArray[s->atlasOff[p]];
2860: if (!cDim) {
2861: if (orientation >= 0) {
2862: const PetscInt dim = s->atlasDof[p];
2863: PetscInt i;
2865: if (mode == INSERT_VALUES) {
2866: for (i = 0; i < dim; ++i) array[i] = values ? values[i] : i;
2867: } else {
2868: for (i = 0; i < dim; ++i) array[i] += values[i];
2869: }
2870: } else {
2871: PetscInt offset = 0;
2872: PetscInt j = -1, field, i;
2874: for (field = 0; field < s->numFields; ++field) {
2875: const PetscInt dim = s->field[field]->atlasDof[p];
2877: for (i = dim - 1; i >= 0; --i) array[++j] = values ? values[i + offset] : i + offset;
2878: offset += dim;
2879: }
2880: }
2881: } else {
2882: if (orientation >= 0) {
2883: const PetscInt dim = s->atlasDof[p];
2884: PetscInt cInd = 0, i;
2885: const PetscInt *cDof;
2887: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2888: if (mode == INSERT_VALUES) {
2889: for (i = 0; i < dim; ++i) {
2890: if ((cInd < cDim) && (i == cDof[cInd])) {
2891: ++cInd;
2892: continue;
2893: }
2894: array[i] = values ? values[i] : i;
2895: }
2896: } else {
2897: for (i = 0; i < dim; ++i) {
2898: if ((cInd < cDim) && (i == cDof[cInd])) {
2899: ++cInd;
2900: continue;
2901: }
2902: array[i] += values[i];
2903: }
2904: }
2905: } else {
2906: const PetscInt *cDof;
2907: PetscInt offset = 0;
2908: PetscInt cOffset = 0;
2909: PetscInt j = 0, field;
2911: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2912: for (field = 0; field < s->numFields; ++field) {
2913: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2914: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2915: const PetscInt sDim = dim - tDim;
2916: PetscInt cInd = 0, i, k;
2918: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2919: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2920: ++cInd;
2921: continue;
2922: }
2923: array[j] = values ? values[k] : k;
2924: }
2925: offset += dim;
2926: cOffset += dim - tDim;
2927: }
2928: }
2929: }
2930: PetscFunctionReturn(PETSC_SUCCESS);
2931: }
2933: /*@
2934: PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs
2936: Not Collective
2938: Input Parameter:
2939: . s - The `PetscSection`
2941: Output Parameter:
2942: . hasConstraints - flag indicating that the section has constrained dofs
2944: Level: intermediate
2946: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2947: @*/
2948: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2949: {
2950: PetscFunctionBegin;
2952: PetscAssertPointer(hasConstraints, 2);
2953: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2954: PetscFunctionReturn(PETSC_SUCCESS);
2955: }
2957: /*@C
2958: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained for a given point
2960: Not Collective
2962: Input Parameters:
2963: + s - The `PetscSection`
2964: - point - The point
2966: Output Parameter:
2967: . indices - The constrained dofs
2969: Level: intermediate
2971: Fortran Notes:
2972: Use `PetscSectionRestoreConstraintIndices()` when the indices are no longer needed
2974: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2975: @*/
2976: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt *indices[])
2977: {
2978: PetscFunctionBegin;
2980: if (s->bc) PetscCall(VecIntGetValuesSection_Private(s->bcIndices, s->bc, point, indices));
2981: else *indices = NULL;
2982: PetscFunctionReturn(PETSC_SUCCESS);
2983: }
2985: /*@
2986: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2988: Not Collective
2990: Input Parameters:
2991: + s - The `PetscSection`
2992: . point - The point
2993: - indices - The constrained dofs
2995: Level: intermediate
2997: .seealso: [PetscSection](ch_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2998: @*/
2999: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
3000: {
3001: PetscFunctionBegin;
3003: if (s->bc) {
3004: const PetscInt dof = s->atlasDof[point];
3005: const PetscInt cdof = s->bc->atlasDof[point];
3006: PetscInt d;
3008: if (indices)
3009: 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]);
3010: PetscCall(VecIntSetValuesSection_Private(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
3011: }
3012: PetscFunctionReturn(PETSC_SUCCESS);
3013: }
3015: /*@C
3016: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
3018: Not Collective
3020: Input Parameters:
3021: + s - The `PetscSection`
3022: . field - The field number
3023: - point - The point
3025: Output Parameter:
3026: . indices - The constrained dofs sorted in ascending order, the length is returned by `PetscSectionGetConstraintDof()`.
3028: Level: intermediate
3030: Fortran Notes:
3031: Use `PetscSectionRestoreFieldConstraintIndices()` to restore the indices when no longer needed
3033: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
3034: @*/
3035: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt *indices[])
3036: {
3037: PetscFunctionBegin;
3039: PetscAssertPointer(indices, 4);
3040: PetscSectionCheckValidField(field, s->numFields);
3041: PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
3042: PetscFunctionReturn(PETSC_SUCCESS);
3043: }
3045: /*@
3046: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
3048: Not Collective
3050: Input Parameters:
3051: + s - The `PetscSection`
3052: . point - The point
3053: . field - The field number
3054: - indices - The constrained dofs
3056: Level: intermediate
3058: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
3059: @*/
3060: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
3061: {
3062: PetscFunctionBegin;
3064: PetscSectionCheckValidField(field, s->numFields);
3065: PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
3066: PetscFunctionReturn(PETSC_SUCCESS);
3067: }
3069: /*@
3070: PetscSectionPermute - Reorder the section according to the input point permutation
3072: Collective
3074: Input Parameters:
3075: + section - The `PetscSection` object
3076: - permutation - The point permutation, old point p becomes new point perm[p]
3078: Output Parameter:
3079: . sectionNew - The permuted `PetscSection`
3081: Level: intermediate
3083: Note:
3084: The data and the access to the data via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` are both changed in `sectionNew`
3086: Compare to `PetscSectionSetPermutation()`
3088: .seealso: [PetscSection](ch_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
3089: @*/
3090: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
3091: {
3092: PetscSection s = section, sNew;
3093: const PetscInt *perm;
3094: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
3096: PetscFunctionBegin;
3099: PetscAssertPointer(sectionNew, 3);
3100: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
3101: PetscCall(PetscSectionGetNumFields(s, &numFields));
3102: if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
3103: for (f = 0; f < numFields; ++f) {
3104: const char *name;
3105: PetscInt numComp;
3107: PetscCall(PetscSectionGetFieldName(s, f, &name));
3108: PetscCall(PetscSectionSetFieldName(sNew, f, name));
3109: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
3110: PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
3111: for (c = 0; c < s->numFieldComponents[f]; ++c) {
3112: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
3113: PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
3114: }
3115: }
3116: PetscCall(ISGetLocalSize(permutation, &numPoints));
3117: PetscCall(ISGetIndices(permutation, &perm));
3118: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
3119: PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
3120: PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
3121: for (p = pStart; p < pEnd; ++p) {
3122: PetscInt dof, cdof;
3124: PetscCall(PetscSectionGetDof(s, p, &dof));
3125: PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
3126: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
3127: if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
3128: for (f = 0; f < numFields; ++f) {
3129: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
3130: PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
3131: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
3132: if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
3133: }
3134: }
3135: PetscCall(PetscSectionSetUp(sNew));
3136: for (p = pStart; p < pEnd; ++p) {
3137: const PetscInt *cind;
3138: PetscInt cdof;
3140: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
3141: if (cdof) {
3142: PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
3143: PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
3144: }
3145: for (f = 0; f < numFields; ++f) {
3146: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
3147: if (cdof) {
3148: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
3149: PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
3150: }
3151: }
3152: }
3153: PetscCall(ISRestoreIndices(permutation, &perm));
3154: *sectionNew = sNew;
3155: PetscFunctionReturn(PETSC_SUCCESS);
3156: }
3158: /*@
3159: PetscSectionSetClosureIndex - Create an internal data structure to speed up closure queries.
3161: Collective
3163: Input Parameters:
3164: + section - The `PetscSection`
3165: . obj - A `PetscObject` which serves as the key for this index
3166: . clSection - `PetscSection` giving the size of the closure of each point
3167: - clPoints - `IS` giving the points in each closure
3169: Level: advanced
3171: Note:
3172: This function creates an internal map from each point to its closure. We compress out closure points with no dofs in this section.
3174: Developer Notes:
3175: The information provided here is completely opaque
3177: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
3178: @*/
3179: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
3180: {
3181: PetscFunctionBegin;
3185: if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
3186: section->clObj = obj;
3187: PetscCall(PetscObjectReference((PetscObject)clSection));
3188: PetscCall(PetscObjectReference((PetscObject)clPoints));
3189: PetscCall(PetscSectionDestroy(§ion->clSection));
3190: PetscCall(ISDestroy(§ion->clPoints));
3191: section->clSection = clSection;
3192: section->clPoints = clPoints;
3193: PetscFunctionReturn(PETSC_SUCCESS);
3194: }
3196: /*@
3197: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
3199: Collective
3201: Input Parameters:
3202: + section - The `PetscSection`
3203: - obj - A `PetscObject` which serves as the key for this index
3205: Output Parameters:
3206: + clSection - `PetscSection` giving the size of the closure of each point
3207: - clPoints - `IS` giving the points in each closure
3209: Level: advanced
3211: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3212: @*/
3213: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
3214: {
3215: PetscFunctionBegin;
3216: if (section->clObj == obj) {
3217: if (clSection) *clSection = section->clSection;
3218: if (clPoints) *clPoints = section->clPoints;
3219: } else {
3220: if (clSection) *clSection = NULL;
3221: if (clPoints) *clPoints = NULL;
3222: }
3223: PetscFunctionReturn(PETSC_SUCCESS);
3224: }
3226: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
3227: {
3228: PetscInt i;
3229: khiter_t iter;
3230: int new_entry;
3231: PetscSectionClosurePermKey key = {depth, clSize};
3232: PetscSectionClosurePermVal *val;
3234: PetscFunctionBegin;
3235: if (section->clObj != obj) {
3236: PetscCall(PetscSectionDestroy(§ion->clSection));
3237: PetscCall(ISDestroy(§ion->clPoints));
3238: }
3239: section->clObj = obj;
3240: if (!section->clHash) PetscCall(PetscClPermCreate(§ion->clHash));
3241: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
3242: val = &kh_val(section->clHash, iter);
3243: if (!new_entry) {
3244: PetscCall(PetscFree(val->perm));
3245: PetscCall(PetscFree(val->invPerm));
3246: }
3247: if (mode == PETSC_COPY_VALUES) {
3248: PetscCall(PetscMalloc1(clSize, &val->perm));
3249: PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
3250: } else if (mode == PETSC_OWN_POINTER) {
3251: val->perm = clPerm;
3252: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
3253: PetscCall(PetscMalloc1(clSize, &val->invPerm));
3254: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
3255: PetscFunctionReturn(PETSC_SUCCESS);
3256: }
3258: /*@
3259: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3261: Not Collective
3263: Input Parameters:
3264: + section - The `PetscSection`
3265: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3266: . depth - Depth of points on which to apply the given permutation
3267: - perm - Permutation of the cell dof closure
3269: Level: intermediate
3271: Notes:
3272: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
3273: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
3274: each topology and degree.
3276: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
3277: exotic/enriched spaces on mixed topology meshes.
3279: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
3280: @*/
3281: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
3282: {
3283: const PetscInt *clPerm = NULL;
3284: PetscInt clSize = 0;
3286: PetscFunctionBegin;
3287: if (perm) {
3288: PetscCall(ISGetLocalSize(perm, &clSize));
3289: PetscCall(ISGetIndices(perm, &clPerm));
3290: }
3291: PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
3292: if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
3293: PetscFunctionReturn(PETSC_SUCCESS);
3294: }
3296: static PetscErrorCode PetscSectionGetClosurePermutation_Private(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3297: {
3298: PetscFunctionBegin;
3299: if (section->clObj == obj) {
3300: PetscSectionClosurePermKey k = {depth, size};
3301: PetscSectionClosurePermVal v;
3303: PetscCall(PetscClPermGet(section->clHash, k, &v));
3304: if (perm) *perm = v.perm;
3305: } else {
3306: if (perm) *perm = NULL;
3307: }
3308: PetscFunctionReturn(PETSC_SUCCESS);
3309: }
3311: /*@
3312: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3314: Not Collective
3316: Input Parameters:
3317: + section - The `PetscSection`
3318: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
3319: . depth - Depth stratum on which to obtain closure permutation
3320: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3322: Output Parameter:
3323: . perm - The dof closure permutation
3325: Level: intermediate
3327: Note:
3328: The user must destroy the `IS` that is returned.
3330: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3331: @*/
3332: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3333: {
3334: const PetscInt *clPerm = NULL;
3336: PetscFunctionBegin;
3337: PetscCall(PetscSectionGetClosurePermutation_Private(section, obj, depth, clSize, &clPerm));
3338: PetscCheck(clPerm, PetscObjectComm(obj), PETSC_ERR_ARG_WRONG, "There is no closure permutation associated with this object for depth %" PetscInt_FMT " of size %" PetscInt_FMT, depth, clSize);
3339: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3340: PetscFunctionReturn(PETSC_SUCCESS);
3341: }
3343: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3344: {
3345: PetscFunctionBegin;
3346: if (section->clObj == obj && section->clHash) {
3347: PetscSectionClosurePermKey k = {depth, size};
3348: PetscSectionClosurePermVal v;
3349: PetscCall(PetscClPermGet(section->clHash, k, &v));
3350: if (perm) *perm = v.invPerm;
3351: } else {
3352: if (perm) *perm = NULL;
3353: }
3354: PetscFunctionReturn(PETSC_SUCCESS);
3355: }
3357: /*@
3358: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
3360: Not Collective
3362: Input Parameters:
3363: + section - The `PetscSection`
3364: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3365: . depth - Depth stratum on which to obtain closure permutation
3366: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3368: Output Parameter:
3369: . perm - The dof closure permutation
3371: Level: intermediate
3373: Note:
3374: The user must destroy the `IS` that is returned.
3376: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3377: @*/
3378: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3379: {
3380: const PetscInt *clPerm = NULL;
3382: PetscFunctionBegin;
3383: PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3384: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3385: PetscFunctionReturn(PETSC_SUCCESS);
3386: }
3388: /*@
3389: PetscSectionGetField - Get the `PetscSection` associated with a single field
3391: Input Parameters:
3392: + s - The `PetscSection`
3393: - field - The field number
3395: Output Parameter:
3396: . subs - The `PetscSection` for the given field, note the chart of `subs` is not set
3398: Level: intermediate
3400: Note:
3401: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
3403: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3404: @*/
3405: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3406: {
3407: PetscFunctionBegin;
3409: PetscAssertPointer(subs, 3);
3410: PetscSectionCheckValidField(field, s->numFields);
3411: *subs = s->field[field];
3412: PetscFunctionReturn(PETSC_SUCCESS);
3413: }
3415: PetscClassId PETSC_SECTION_SYM_CLASSID;
3416: PetscFunctionList PetscSectionSymList = NULL;
3418: /*@
3419: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
3421: Collective
3423: Input Parameter:
3424: . comm - the MPI communicator
3426: Output Parameter:
3427: . sym - pointer to the new set of symmetries
3429: Level: developer
3431: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3432: @*/
3433: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3434: {
3435: PetscFunctionBegin;
3436: PetscAssertPointer(sym, 2);
3437: PetscCall(ISInitializePackage());
3439: PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3440: PetscFunctionReturn(PETSC_SUCCESS);
3441: }
3443: /*@
3444: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
3446: Collective
3448: Input Parameters:
3449: + sym - The section symmetry object
3450: - method - The name of the section symmetry type
3452: Level: developer
3454: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3455: @*/
3456: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3457: {
3458: PetscErrorCode (*r)(PetscSectionSym);
3459: PetscBool match;
3461: PetscFunctionBegin;
3463: PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3464: if (match) PetscFunctionReturn(PETSC_SUCCESS);
3466: PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3467: PetscCheck(r, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3468: PetscTryTypeMethod(sym, destroy);
3469: sym->ops->destroy = NULL;
3471: PetscCall((*r)(sym));
3472: PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3473: PetscFunctionReturn(PETSC_SUCCESS);
3474: }
3476: /*@
3477: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
3479: Not Collective
3481: Input Parameter:
3482: . sym - The section symmetry
3484: Output Parameter:
3485: . type - The index set type name
3487: Level: developer
3489: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3490: @*/
3491: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3492: {
3493: PetscFunctionBegin;
3495: PetscAssertPointer(type, 2);
3496: *type = ((PetscObject)sym)->type_name;
3497: PetscFunctionReturn(PETSC_SUCCESS);
3498: }
3500: /*@C
3501: PetscSectionSymRegister - Registers a new section symmetry implementation
3503: Not Collective, No Fortran Support
3505: Input Parameters:
3506: + sname - The name of a new user-defined creation routine
3507: - function - The creation routine itself
3509: Level: developer
3511: Notes:
3512: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
3514: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3515: @*/
3516: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3517: {
3518: PetscFunctionBegin;
3519: PetscCall(ISInitializePackage());
3520: PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3521: PetscFunctionReturn(PETSC_SUCCESS);
3522: }
3524: /*@
3525: PetscSectionSymDestroy - Destroys a section symmetry.
3527: Collective
3529: Input Parameter:
3530: . sym - the section symmetry
3532: Level: developer
3534: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`
3535: @*/
3536: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3537: {
3538: SymWorkLink link, next;
3540: PetscFunctionBegin;
3541: if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3543: if (--((PetscObject)*sym)->refct > 0) {
3544: *sym = NULL;
3545: PetscFunctionReturn(PETSC_SUCCESS);
3546: }
3547: PetscTryTypeMethod(*sym, destroy);
3548: PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3549: for (link = (*sym)->workin; link; link = next) {
3550: PetscInt **perms = (PetscInt **)link->perms;
3551: PetscScalar **rots = (PetscScalar **)link->rots;
3552: PetscCall(PetscFree2(perms, rots));
3553: next = link->next;
3554: PetscCall(PetscFree(link));
3555: }
3556: (*sym)->workin = NULL;
3557: PetscCall(PetscHeaderDestroy(sym));
3558: PetscFunctionReturn(PETSC_SUCCESS);
3559: }
3561: /*@
3562: PetscSectionSymView - Displays a section symmetry
3564: Collective
3566: Input Parameters:
3567: + sym - the index set
3568: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3570: Level: developer
3572: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3573: @*/
3574: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3575: {
3576: PetscFunctionBegin;
3578: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3580: PetscCheckSameComm(sym, 1, viewer, 2);
3581: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3582: PetscTryTypeMethod(sym, view, viewer);
3583: PetscFunctionReturn(PETSC_SUCCESS);
3584: }
3586: /*@
3587: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3589: Collective
3591: Input Parameters:
3592: + section - the section describing data layout
3593: - sym - the symmetry describing the affect of orientation on the access of the data
3595: Level: developer
3597: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3598: @*/
3599: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3600: {
3601: PetscFunctionBegin;
3603: PetscCall(PetscSectionSymDestroy(§ion->sym));
3604: if (sym) {
3606: PetscCheckSameComm(section, 1, sym, 2);
3607: PetscCall(PetscObjectReference((PetscObject)sym));
3608: }
3609: section->sym = sym;
3610: PetscFunctionReturn(PETSC_SUCCESS);
3611: }
3613: /*@
3614: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3616: Not Collective
3618: Input Parameter:
3619: . section - the section describing data layout
3621: Output Parameter:
3622: . sym - the symmetry describing the affect of orientation on the access of the data, provided previously by `PetscSectionSetSym()`
3624: Level: developer
3626: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3627: @*/
3628: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3629: {
3630: PetscFunctionBegin;
3632: *sym = section->sym;
3633: PetscFunctionReturn(PETSC_SUCCESS);
3634: }
3636: /*@
3637: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3639: Collective
3641: Input Parameters:
3642: + section - the section describing data layout
3643: . field - the field number
3644: - sym - the symmetry describing the affect of orientation on the access of the data
3646: Level: developer
3648: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3649: @*/
3650: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3651: {
3652: PetscFunctionBegin;
3654: PetscSectionCheckValidField(field, section->numFields);
3655: PetscCall(PetscSectionSetSym(section->field[field], sym));
3656: PetscFunctionReturn(PETSC_SUCCESS);
3657: }
3659: /*@
3660: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3662: Collective
3664: Input Parameters:
3665: + section - the section describing data layout
3666: - field - the field number
3668: Output Parameter:
3669: . sym - the symmetry describing the affect of orientation on the access of the data
3671: Level: developer
3673: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3674: @*/
3675: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3676: {
3677: PetscFunctionBegin;
3679: PetscSectionCheckValidField(field, section->numFields);
3680: *sym = section->field[field]->sym;
3681: PetscFunctionReturn(PETSC_SUCCESS);
3682: }
3684: /*@C
3685: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3687: Not Collective
3689: Input Parameters:
3690: + section - the section
3691: . numPoints - the number of points
3692: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3693: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3694: context, see `DMPlexGetConeOrientation()`).
3696: Output Parameters:
3697: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3698: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3699: identity).
3701: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3702: .vb
3703: const PetscInt **perms;
3704: const PetscScalar **rots;
3705: PetscInt lOffset;
3707: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3708: for (i = 0, lOffset = 0; i < numPoints; i++) {
3709: PetscInt point = points[2*i], dof, sOffset;
3710: const PetscInt *perm = perms ? perms[i] : NULL;
3711: const PetscScalar *rot = rots ? rots[i] : NULL;
3713: PetscSectionGetDof(section,point,&dof);
3714: PetscSectionGetOffset(section,point,&sOffset);
3716: if (perm) { for (j = 0; j < dof; j++) lArray[lOffset + perm[j]] = sArray[sOffset + j]; }
3717: else { for (j = 0; j < dof; j++) lArray[lOffset + j ] = sArray[sOffset + j]; }
3718: if (rot) { for (j = 0; j < dof; j++) lArray[lOffset + j ] *= rot[j]; }
3719: lOffset += dof;
3720: }
3721: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3722: .ve
3724: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3725: .vb
3726: const PetscInt **perms;
3727: const PetscScalar **rots;
3728: PetscInt lOffset;
3730: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3731: for (i = 0, lOffset = 0; i < numPoints; i++) {
3732: PetscInt point = points[2*i], dof, sOffset;
3733: const PetscInt *perm = perms ? perms[i] : NULL;
3734: const PetscScalar *rot = rots ? rots[i] : NULL;
3736: PetscSectionGetDof(section,point,&dof);
3737: PetscSectionGetOffset(section,point,&sOff);
3739: if (perm) { for (j = 0; j < dof; j++) sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.); }
3740: else { for (j = 0; j < dof; j++) sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.); }
3741: offset += dof;
3742: }
3743: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3744: .ve
3746: Level: developer
3748: Notes:
3749: `PetscSectionSetSym()` must have been previously called to provide the symmetries to the `PetscSection`
3751: Use `PetscSectionRestorePointSyms()` when finished with the data
3753: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3754: @*/
3755: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3756: {
3757: PetscSectionSym sym;
3759: PetscFunctionBegin;
3761: if (numPoints) PetscAssertPointer(points, 3);
3762: if (perms) *perms = NULL;
3763: if (rots) *rots = NULL;
3764: sym = section->sym;
3765: if (sym && (perms || rots)) {
3766: SymWorkLink link;
3768: if (sym->workin) {
3769: link = sym->workin;
3770: sym->workin = sym->workin->next;
3771: } else {
3772: PetscCall(PetscNew(&link));
3773: }
3774: if (numPoints > link->numPoints) {
3775: PetscInt **perms = (PetscInt **)link->perms;
3776: PetscScalar **rots = (PetscScalar **)link->rots;
3777: PetscCall(PetscFree2(perms, rots));
3778: PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3779: link->numPoints = numPoints;
3780: }
3781: link->next = sym->workout;
3782: sym->workout = link;
3783: PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3784: PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3785: PetscUseTypeMethod(sym, getpoints, section, numPoints, points, link->perms, link->rots);
3786: if (perms) *perms = link->perms;
3787: if (rots) *rots = link->rots;
3788: }
3789: PetscFunctionReturn(PETSC_SUCCESS);
3790: }
3792: /*@C
3793: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3795: Not Collective
3797: Input Parameters:
3798: + section - the section
3799: . numPoints - the number of points
3800: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3801: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3802: context, see `DMPlexGetConeOrientation()`).
3803: . perms - The permutations for the given orientations: set to `NULL` at conclusion
3804: - rots - The field rotations symmetries for the given orientations: set to `NULL` at conclusion
3806: Level: developer
3808: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3809: @*/
3810: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3811: {
3812: PetscSectionSym sym;
3814: PetscFunctionBegin;
3816: sym = section->sym;
3817: if (sym && (perms || rots)) {
3818: SymWorkLink *p, link;
3820: for (p = &sym->workout; (link = *p); p = &link->next) {
3821: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3822: *p = link->next;
3823: link->next = sym->workin;
3824: sym->workin = link;
3825: if (perms) *perms = NULL;
3826: if (rots) *rots = NULL;
3827: PetscFunctionReturn(PETSC_SUCCESS);
3828: }
3829: }
3830: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3831: }
3832: PetscFunctionReturn(PETSC_SUCCESS);
3833: }
3835: /*@C
3836: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3838: Not Collective
3840: Input Parameters:
3841: + section - the section
3842: . field - the field of the section
3843: . numPoints - the number of points
3844: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3845: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3846: context, see `DMPlexGetConeOrientation()`).
3848: Output Parameters:
3849: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3850: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3851: identity).
3853: Level: developer
3855: Notes:
3856: `PetscSectionSetFieldSym()` must have been previously called to provide the symmetries to the `PetscSection`
3858: Use `PetscSectionRestoreFieldPointSyms()` when finished with the data
3860: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3861: @*/
3862: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3863: {
3864: PetscFunctionBegin;
3866: 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);
3867: PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3868: PetscFunctionReturn(PETSC_SUCCESS);
3869: }
3871: /*@C
3872: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3874: Not Collective
3876: Input Parameters:
3877: + section - the section
3878: . field - the field number
3879: . numPoints - the number of points
3880: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3881: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3882: context, see `DMPlexGetConeOrientation()`).
3883: . perms - The permutations for the given orientations: set to NULL at conclusion
3884: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3886: Level: developer
3888: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3889: @*/
3890: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3891: {
3892: PetscFunctionBegin;
3894: 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);
3895: PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3896: PetscFunctionReturn(PETSC_SUCCESS);
3897: }
3899: /*@
3900: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3902: Not Collective
3904: Input Parameter:
3905: . sym - the `PetscSectionSym`
3907: Output Parameter:
3908: . nsym - the equivalent symmetries
3910: Level: developer
3912: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3913: @*/
3914: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3915: {
3916: PetscFunctionBegin;
3919: PetscTryTypeMethod(sym, copy, nsym);
3920: PetscFunctionReturn(PETSC_SUCCESS);
3921: }
3923: /*@
3924: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3926: Collective
3928: Input Parameters:
3929: + sym - the `PetscSectionSym`
3930: - migrationSF - the distribution map from roots to leaves
3932: Output Parameter:
3933: . dsym - the redistributed symmetries
3935: Level: developer
3937: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3938: @*/
3939: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3940: {
3941: PetscFunctionBegin;
3944: PetscAssertPointer(dsym, 3);
3945: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3946: PetscFunctionReturn(PETSC_SUCCESS);
3947: }
3949: /*@
3950: PetscSectionGetUseFieldOffsets - Get the flag indicating if field offsets are used directly in a global section, rather than just the point offset
3952: Not Collective
3954: Input Parameter:
3955: . s - the global `PetscSection`
3957: Output Parameter:
3958: . flg - the flag
3960: Level: developer
3962: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3963: @*/
3964: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3965: {
3966: PetscFunctionBegin;
3968: *flg = s->useFieldOff;
3969: PetscFunctionReturn(PETSC_SUCCESS);
3970: }
3972: /*@
3973: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3975: Not Collective
3977: Input Parameters:
3978: + s - the global `PetscSection`
3979: - flg - the flag
3981: Level: developer
3983: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3984: @*/
3985: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3986: {
3987: PetscFunctionBegin;
3989: s->useFieldOff = flg;
3990: PetscFunctionReturn(PETSC_SUCCESS);
3991: }
3993: #define PetscSectionExpandPoints_Loop(TYPE) \
3994: do { \
3995: PetscInt i, n, o0, o1, size; \
3996: TYPE *a0 = (TYPE *)origArray, *a1; \
3997: PetscCall(PetscSectionGetStorageSize(s, &size)); \
3998: PetscCall(PetscMalloc1(size, &a1)); \
3999: for (i = 0; i < npoints; i++) { \
4000: PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
4001: PetscCall(PetscSectionGetOffset(s, i, &o1)); \
4002: PetscCall(PetscSectionGetDof(s, i, &n)); \
4003: PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n * unitsize)); \
4004: } \
4005: *newArray = (void *)a1; \
4006: } while (0)
4008: /*@
4009: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
4011: Not Collective
4013: Input Parameters:
4014: + origSection - the `PetscSection` describing the layout of the array
4015: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
4016: . origArray - the array; its size must be equal to the storage size of `origSection`
4017: - points - `IS` with points to extract; its indices must lie in the chart of `origSection`
4019: Output Parameters:
4020: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
4021: - newArray - the array of the extracted DOFs; its size is the storage size of `newSection`
4023: Level: developer
4025: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
4026: @*/
4027: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
4028: {
4029: PetscSection s;
4030: const PetscInt *points_;
4031: PetscInt i, n, npoints, pStart, pEnd;
4032: PetscMPIInt unitsize;
4034: PetscFunctionBegin;
4036: PetscAssertPointer(origArray, 3);
4038: if (newSection) PetscAssertPointer(newSection, 5);
4039: if (newArray) PetscAssertPointer(newArray, 6);
4040: PetscCallMPI(MPI_Type_size(dataType, &unitsize));
4041: PetscCall(ISGetLocalSize(points, &npoints));
4042: PetscCall(ISGetIndices(points, &points_));
4043: PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
4044: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
4045: PetscCall(PetscSectionSetChart(s, 0, npoints));
4046: for (i = 0; i < npoints; i++) {
4047: 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);
4048: PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
4049: PetscCall(PetscSectionSetDof(s, i, n));
4050: }
4051: PetscCall(PetscSectionSetUp(s));
4052: if (newArray) {
4053: if (dataType == MPIU_INT) {
4054: PetscSectionExpandPoints_Loop(PetscInt);
4055: } else if (dataType == MPIU_SCALAR) {
4056: PetscSectionExpandPoints_Loop(PetscScalar);
4057: } else if (dataType == MPIU_REAL) {
4058: PetscSectionExpandPoints_Loop(PetscReal);
4059: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
4060: }
4061: if (newSection) {
4062: *newSection = s;
4063: } else {
4064: PetscCall(PetscSectionDestroy(&s));
4065: }
4066: PetscCall(ISRestoreIndices(points, &points_));
4067: PetscFunctionReturn(PETSC_SUCCESS);
4068: }
4070: /*@
4071: PetscSectionMigrateData - Migrate data described by a `PetscSection` using a `PetscSF` that defines a original-to-new (root-to-leaf) point mapping
4073: Collective
4075: Input Parameters:
4076: + migratePointSF - defines the mapping (communication) of the root points to the leaf points
4077: . datatype - the type of data
4078: . rootSection - the `PetscSection` that describes the data layout on the root points (how many dof and what fields are associated with each root point)
4079: - rootData - the existing data array described by `rootSection`, may be `NULL` is storage size of `rootSection` is zero
4081: Output Parameters:
4082: + leafSection - the new `PetscSection` that describes the data layout on the leaf points
4083: . leafData - the redistributed data array that is associated with the leaf points
4084: - migrateDataSF - defines the mapping (communication) of the `rootData` array to the `leafData` array, may be `NULL` if not needed
4086: Level: advanced
4088: Notes:
4089: This function can best be thought of as applying `PetscSFBcastBegin()` to an array described by a `PetscSection`.
4090: While `PetscSFBcastBegin()` is limited to broadcasting data that is of the same size for every index, this function allows the data to be a different size for each index.
4091: The size and layout of that irregularly sized data before and after `PetscSFBcastBegin()` is described by the `rootSection` and `leafSection`, respectively.
4093: This function combines `PetscSFDistributeSection()`, `PetscSFCreateSectionSF()`, and `PetscSFBcastBegin()`/`PetscSFBcastEnd()` into a single call.
4094: `migrateDataSF` can be used to repeat the `PetscSFBcastBegin()`/`PetscSFBcastEnd()` on a different data array described by the same `rootSection`.
4096: This should not be used for global-to-local type communciation patterns.
4097: For this use case, see `PetscSectionCreateGlobalSection()` and `PetscSFSetGraphSection()`.
4099: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSFDistributeSection()`, `PetscSFCreateSectionSF()`, `DMPlexDistributeData()`
4100: @*/
4101: PetscErrorCode PetscSectionMigrateData(PetscSF migratePointSF, MPI_Datatype datatype, PetscSection rootSection, const void *rootData, PetscSection leafSection, void *leafData[], PetscSF *migrateDataSF)
4102: {
4103: PetscSF fieldSF;
4104: PetscInt *remoteOffsets, fieldSize;
4105: PetscMPIInt dataSize;
4107: PetscFunctionBegin;
4110: if (rootData) PetscAssertPointer(rootData, 4);
4111: else {
4112: PetscInt size;
4113: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
4114: PetscCheck(size == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "originalData may be NULL iff the storage size of originalSection is zero, but is %" PetscInt_FMT, size);
4115: }
4117: PetscAssertPointer(leafData, 6);
4118: if (migrateDataSF) PetscAssertPointer(migrateDataSF, 7);
4120: PetscCall(PetscSFDistributeSection(migratePointSF, rootSection, &remoteOffsets, leafSection));
4121: PetscCall(PetscSFCreateSectionSF(migratePointSF, rootSection, remoteOffsets, leafSection, &fieldSF));
4122: PetscCall(PetscFree(remoteOffsets));
4124: PetscCall(PetscSectionGetStorageSize(leafSection, &fieldSize));
4125: PetscCallMPI(MPI_Type_size(datatype, &dataSize));
4126: PetscCall(PetscMalloc(fieldSize * dataSize, leafData));
4127: PetscCall(PetscSFBcastBegin(fieldSF, datatype, rootData, *leafData, MPI_REPLACE));
4128: PetscCall(PetscSFBcastEnd(fieldSF, datatype, rootData, *leafData, MPI_REPLACE));
4130: if (migrateDataSF) *migrateDataSF = fieldSF;
4131: else PetscCall(PetscSFDestroy(&fieldSF));
4132: PetscFunctionReturn(PETSC_SUCCESS);
4133: }