Actual source code: section.c
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsf.h>
8: PetscClassId PETSC_SECTION_CLASSID;
10: /*@
11: PetscSectionCreate - Allocates a `PetscSection` and sets the map contents to the default.
13: Collective
15: Input Parameters:
16: + comm - the MPI communicator
17: - s - pointer to the section
19: Level: beginner
21: Notes:
22: Typical calling sequence
23: .vb
24: PetscSectionCreate(MPI_Comm,PetscSection *);!
25: PetscSectionSetNumFields(PetscSection, numFields);
26: PetscSectionSetChart(PetscSection,low,high);
27: PetscSectionSetDof(PetscSection,point,numdof);
28: PetscSectionSetUp(PetscSection);
29: PetscSectionGetOffset(PetscSection,point,PetscInt *);
30: PetscSectionDestroy(PetscSection);
31: .ve
33: The `PetscSection` object and methods are intended to be used in the PETSc `Vec` and `Mat` implementations. The indices returned by the `PetscSection` are appropriate for the kind of `Vec` it is associated with. For example, if the vector being indexed is a local vector, we call the section a local section. If the section indexes a global vector, we call it a global section. For parallel vectors, like global vectors, we use negative indices to indicate dofs owned by other processes.
35: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionDestroy()`, `PetscSectionCreateGlobalSection()`
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
39: PetscFunctionBegin;
40: PetscAssertPointer(s, 2);
41: PetscCall(ISInitializePackage());
43: PetscCall(PetscHeaderCreate(*s, PETSC_SECTION_CLASSID, "PetscSection", "Section", "IS", comm, PetscSectionDestroy, PetscSectionView));
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](sec_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](sec_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](sec_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](sec_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, MPIU_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](sec_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](sec_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, (char **)&s->fieldNames[f]));
399: PetscCall(PetscSNPrintf(name, 64, "Component_0"));
400: PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]));
401: PetscCall(PetscStrallocpy(name, (char **)&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](sec_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](sec_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, (char **)&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](sec_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](sec_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, (char **)&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](sec_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](sec_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, (char ***)&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, (char **)&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: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionCreate()`
628: @*/
629: PetscErrorCode PetscSectionGetChart(PetscSection s, PetscInt *pStart, PetscInt *pEnd)
630: {
631: PetscFunctionBegin;
633: if (pStart) *pStart = s->pStart;
634: if (pEnd) *pEnd = s->pEnd;
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: /*@
639: PetscSectionSetChart - Sets the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
641: Not Collective
643: Input Parameters:
644: + s - the `PetscSection`
645: . pStart - the first point
646: - pEnd - one past the last point, `pStart` $ \le $ `pEnd`
648: Level: intermediate
650: Notes:
651: The charts on different MPI processes may (and often do) overlap
653: If you intend to use `PetscSectionSetNumFields()` it must be called before this call.
655: The chart for all fields created with `PetscSectionSetNumFields()` is the same as this chart.
657: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`, `PetscSectionSetNumFields()`
658: @*/
659: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
660: {
661: PetscInt f;
663: PetscFunctionBegin;
665: PetscCheck(pEnd >= pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Chart pEnd %" PetscInt_FMT " cannot be smaller than chart pStart %" PetscInt_FMT, pEnd, pStart);
666: if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(PETSC_SUCCESS);
667: /* Cannot Reset() because it destroys field information */
668: s->setup = PETSC_FALSE;
669: PetscCall(PetscSectionDestroy(&s->bc));
670: PetscCall(PetscFree(s->bcIndices));
671: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
673: s->pStart = pStart;
674: s->pEnd = pEnd;
675: PetscCall(PetscMalloc2(pEnd - pStart, &s->atlasDof, pEnd - pStart, &s->atlasOff));
676: PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart));
677: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: /*@
682: PetscSectionGetPermutation - Returns the permutation of [0, `pEnd` - `pStart`) or `NULL` that was set with `PetscSectionSetPermutation()`
684: Not Collective
686: Input Parameter:
687: . s - the `PetscSection`
689: Output Parameter:
690: . perm - The permutation as an `IS`
692: Level: intermediate
694: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
695: @*/
696: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
697: {
698: PetscFunctionBegin;
700: if (perm) {
701: PetscAssertPointer(perm, 2);
702: *perm = s->perm;
703: }
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: /*@
708: PetscSectionSetPermutation - Sets a permutation of the chart for this section, [0, `pEnd` - `pStart`), which determines the order to store the `PetscSection` information
710: Not Collective
712: Input Parameters:
713: + s - the `PetscSection`
714: - perm - the permutation of points
716: Level: intermediate
718: Notes:
719: The permutation must be provided before `PetscSectionSetUp()`.
721: The data in the `PetscSection` are permuted but the access via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` is not changed
723: Compare to `PetscSectionPermute()`
725: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetUp()`, `PetscSectionGetPermutation()`, `PetscSectionPermute()`, `PetscSectionCreate()`
726: @*/
727: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
728: {
729: PetscFunctionBegin;
732: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
733: if (s->perm != perm) {
734: PetscCall(ISDestroy(&s->perm));
735: if (perm) {
736: s->perm = perm;
737: PetscCall(PetscObjectReference((PetscObject)s->perm));
738: }
739: }
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@C
744: PetscSectionGetBlockStarts - Returns a table indicating which points start new blocks
746: Not Collective, No Fortran Support
748: Input Parameter:
749: . s - the `PetscSection`
751: Output Parameter:
752: . blockStarts - The `PetscBT` with a 1 for each point that begins a block
754: Notes:
755: The table is on [0, `pEnd` - `pStart`).
757: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
759: Level: intermediate
761: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
762: @*/
763: PetscErrorCode PetscSectionGetBlockStarts(PetscSection s, PetscBT *blockStarts)
764: {
765: PetscFunctionBegin;
767: if (blockStarts) {
768: PetscAssertPointer(blockStarts, 2);
769: *blockStarts = s->blockStarts;
770: }
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: /*@C
775: PetscSectionSetBlockStarts - Sets a table indicating which points start new blocks
777: Not Collective, No Fortran Support
779: Input Parameters:
780: + s - the `PetscSection`
781: - blockStarts - The `PetscBT` with a 1 for each point that begins a block
783: Level: intermediate
785: Notes:
786: The table is on [0, `pEnd` - `pStart`). PETSc takes ownership of the `PetscBT` when it is passed in and will destroy it. The user should not destroy it.
788: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
790: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionGetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
791: @*/
792: PetscErrorCode PetscSectionSetBlockStarts(PetscSection s, PetscBT blockStarts)
793: {
794: PetscFunctionBegin;
796: if (s->blockStarts != blockStarts) {
797: PetscCall(PetscBTDestroy(&s->blockStarts));
798: s->blockStarts = blockStarts;
799: }
800: PetscFunctionReturn(PETSC_SUCCESS);
801: }
803: /*@
804: PetscSectionGetPointMajor - Returns the flag for dof ordering, `PETSC_TRUE` if it is point major, `PETSC_FALSE` if it is field major
806: Not Collective
808: Input Parameter:
809: . s - the `PetscSection`
811: Output Parameter:
812: . pm - the flag for point major ordering
814: Level: intermediate
816: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetPointMajor()`
817: @*/
818: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
819: {
820: PetscFunctionBegin;
822: PetscAssertPointer(pm, 2);
823: *pm = s->pointMajor;
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: /*@
828: PetscSectionSetPointMajor - Sets the flag for dof ordering, `PETSC_TRUE` for point major, otherwise it will be field major
830: Not Collective
832: Input Parameters:
833: + s - the `PetscSection`
834: - pm - the flag for point major ordering
836: Level: intermediate
838: Note:
839: Field-major order is not recommended unless you are managing the entire problem yourself, since many higher-level functions in PETSc depend on point-major order.
841: Point major order means the degrees of freedom are stored as follows
842: .vb
843: all the degrees of freedom for each point are stored contiguously, one point after another (respecting a permutation set with `PetscSectionSetPermutation()`)
844: for each point
845: the degrees of freedom for each field (starting with the unnamed default field) are listed in order by field
846: .ve
848: Field major order means the degrees of freedom are stored as follows
849: .vb
850: all degrees of freedom for each field (including the unnamed default field) are stored contiguously, one field after another
851: for each field (started with unnamed default field)
852: the degrees of freedom for each point are listed in order by point (respecting a permutation set with `PetscSectionSetPermutation()`)
853: .ve
855: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`, `PetscSectionSetPermutation()`
856: @*/
857: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
858: {
859: PetscFunctionBegin;
861: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
862: s->pointMajor = pm;
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: /*@
867: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets in the `PetscSection`.
868: The value is set with `PetscSectionSetIncludesConstraints()`
870: Not Collective
872: Input Parameter:
873: . s - the `PetscSection`
875: Output Parameter:
876: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
878: Level: intermediate
880: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
881: @*/
882: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
883: {
884: PetscFunctionBegin;
886: PetscAssertPointer(includesConstraints, 2);
887: *includesConstraints = s->includesConstraints;
888: PetscFunctionReturn(PETSC_SUCCESS);
889: }
891: /*@
892: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
894: Not Collective
896: Input Parameters:
897: + s - the `PetscSection`
898: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
900: Level: intermediate
902: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
903: @*/
904: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
905: {
906: PetscFunctionBegin;
908: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
909: s->includesConstraints = includesConstraints;
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: /*@
914: PetscSectionGetDof - Return the total number of degrees of freedom associated with a given point.
916: Not Collective
918: Input Parameters:
919: + s - the `PetscSection`
920: - point - the point
922: Output Parameter:
923: . numDof - the number of dof
925: Level: intermediate
927: Notes:
928: In a global section, this size will be negative for points not owned by this process.
930: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
932: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
933: @*/
934: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
935: {
936: PetscFunctionBeginHot;
938: PetscAssertPointer(numDof, 3);
939: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
940: *numDof = s->atlasDof[point - s->pStart];
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: /*@
945: PetscSectionSetDof - Sets the total number of degrees of freedom associated with a given point.
947: Not Collective
949: Input Parameters:
950: + s - the `PetscSection`
951: . point - the point
952: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
954: Level: intermediate
956: Note:
957: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
959: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
960: @*/
961: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
962: {
963: PetscFunctionBegin;
965: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
966: s->atlasDof[point - s->pStart] = numDof;
967: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@
972: PetscSectionAddDof - Adds to the total number of degrees of freedom associated with a given point.
974: Not Collective
976: Input Parameters:
977: + s - the `PetscSection`
978: . point - the point
979: - numDof - the number of additional dof
981: Level: intermediate
983: Note:
984: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
986: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
987: @*/
988: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
989: {
990: PetscFunctionBeginHot;
992: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
993: PetscCheck(numDof >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numDof %" PetscInt_FMT " should not be negative", numDof);
994: s->atlasDof[point - s->pStart] += numDof;
995: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: /*@
1000: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
1002: Not Collective
1004: Input Parameters:
1005: + s - the `PetscSection`
1006: . point - the point
1007: - field - the field
1009: Output Parameter:
1010: . numDof - the number of dof
1012: Level: intermediate
1014: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
1015: @*/
1016: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1017: {
1018: PetscFunctionBegin;
1020: PetscAssertPointer(numDof, 4);
1021: PetscSectionCheckValidField(field, s->numFields);
1022: PetscCall(PetscSectionGetDof(s->field[field], point, numDof));
1023: PetscFunctionReturn(PETSC_SUCCESS);
1024: }
1026: /*@
1027: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
1029: Not Collective
1031: Input Parameters:
1032: + s - the `PetscSection`
1033: . point - the point
1034: . field - the field
1035: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
1037: Level: intermediate
1039: Note:
1040: When setting the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
1041: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1043: This is equivalent to
1044: .vb
1045: PetscSection fs;
1046: PetscSectionGetField(s,field,&fs)
1047: PetscSectionSetDof(fs,numDof)
1048: .ve
1050: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`, `PetscSectionAddDof()`, `PetscSectionSetDof()`
1051: @*/
1052: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1053: {
1054: PetscFunctionBegin;
1056: PetscSectionCheckValidField(field, s->numFields);
1057: PetscCall(PetscSectionSetDof(s->field[field], point, numDof));
1058: PetscFunctionReturn(PETSC_SUCCESS);
1059: }
1061: /*@
1062: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
1064: Not Collective
1066: Input Parameters:
1067: + s - the `PetscSection`
1068: . point - the point
1069: . field - the field
1070: - numDof - the number of dof
1072: Level: intermediate
1074: Notes:
1075: When adding to the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
1076: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1078: This is equivalent to
1079: .vb
1080: PetscSection fs;
1081: PetscSectionGetField(s,field,&fs)
1082: PetscSectionAddDof(fs,numDof)
1083: .ve
1085: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
1086: @*/
1087: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1088: {
1089: PetscFunctionBegin;
1091: PetscSectionCheckValidField(field, s->numFields);
1092: PetscCall(PetscSectionAddDof(s->field[field], point, numDof));
1093: PetscFunctionReturn(PETSC_SUCCESS);
1094: }
1096: /*@
1097: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
1099: Not Collective
1101: Input Parameters:
1102: + s - the `PetscSection`
1103: - point - the point
1105: Output Parameter:
1106: . numDof - the number of dof which are fixed by constraints
1108: Level: intermediate
1110: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
1111: @*/
1112: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
1113: {
1114: PetscFunctionBegin;
1116: PetscAssertPointer(numDof, 3);
1117: if (s->bc) {
1118: PetscCall(PetscSectionGetDof(s->bc, point, numDof));
1119: } else *numDof = 0;
1120: PetscFunctionReturn(PETSC_SUCCESS);
1121: }
1123: /*@
1124: PetscSectionSetConstraintDof - Set the number of constrained degrees of freedom associated with a given point.
1126: Not Collective
1128: Input Parameters:
1129: + s - the `PetscSection`
1130: . point - the point
1131: - numDof - the number of dof which are fixed by constraints
1133: Level: intermediate
1135: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1136: @*/
1137: PetscErrorCode PetscSectionSetConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1138: {
1139: PetscFunctionBegin;
1141: if (numDof) {
1142: PetscCall(PetscSectionCheckConstraints_Private(s));
1143: PetscCall(PetscSectionSetDof(s->bc, point, numDof));
1144: }
1145: PetscFunctionReturn(PETSC_SUCCESS);
1146: }
1148: /*@
1149: PetscSectionAddConstraintDof - Increment the number of constrained degrees of freedom associated with a given point.
1151: Not Collective
1153: Input Parameters:
1154: + s - the `PetscSection`
1155: . point - the point
1156: - numDof - the number of additional dof which are fixed by constraints
1158: Level: intermediate
1160: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetConstraintDof()`, `PetscSectionCreate()`
1161: @*/
1162: PetscErrorCode PetscSectionAddConstraintDof(PetscSection s, PetscInt point, PetscInt numDof)
1163: {
1164: PetscFunctionBegin;
1166: if (numDof) {
1167: PetscCall(PetscSectionCheckConstraints_Private(s));
1168: PetscCall(PetscSectionAddDof(s->bc, point, numDof));
1169: }
1170: PetscFunctionReturn(PETSC_SUCCESS);
1171: }
1173: /*@
1174: PetscSectionGetFieldConstraintDof - Return the number of constrained degrees of freedom associated with a given field on a point.
1176: Not Collective
1178: Input Parameters:
1179: + s - the `PetscSection`
1180: . point - the point
1181: - field - the field
1183: Output Parameter:
1184: . numDof - the number of dof which are fixed by constraints
1186: Level: intermediate
1188: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetFieldConstraintDof()`, `PetscSectionCreate()`
1189: @*/
1190: PetscErrorCode PetscSectionGetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1191: {
1192: PetscFunctionBegin;
1194: PetscAssertPointer(numDof, 4);
1195: PetscSectionCheckValidField(field, s->numFields);
1196: PetscCall(PetscSectionGetConstraintDof(s->field[field], point, numDof));
1197: PetscFunctionReturn(PETSC_SUCCESS);
1198: }
1200: /*@
1201: PetscSectionSetFieldConstraintDof - Set the number of constrained degrees of freedom associated with a given field on a point.
1203: Not Collective
1205: Input Parameters:
1206: + s - the `PetscSection`
1207: . point - the point
1208: . field - the field
1209: - numDof - the number of dof which are fixed by constraints
1211: Level: intermediate
1213: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1214: @*/
1215: PetscErrorCode PetscSectionSetFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1216: {
1217: PetscFunctionBegin;
1219: PetscSectionCheckValidField(field, s->numFields);
1220: PetscCall(PetscSectionSetConstraintDof(s->field[field], point, numDof));
1221: PetscFunctionReturn(PETSC_SUCCESS);
1222: }
1224: /*@
1225: PetscSectionAddFieldConstraintDof - Increment the number of constrained degrees of freedom associated with a given field on a point.
1227: Not Collective
1229: Input Parameters:
1230: + s - the `PetscSection`
1231: . point - the point
1232: . field - the field
1233: - numDof - the number of additional dof which are fixed by constraints
1235: Level: intermediate
1237: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionAddDof()`, `PetscSectionGetFieldConstraintDof()`, `PetscSectionCreate()`
1238: @*/
1239: PetscErrorCode PetscSectionAddFieldConstraintDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1240: {
1241: PetscFunctionBegin;
1243: PetscSectionCheckValidField(field, s->numFields);
1244: PetscCall(PetscSectionAddConstraintDof(s->field[field], point, numDof));
1245: PetscFunctionReturn(PETSC_SUCCESS);
1246: }
1248: /*@
1249: PetscSectionSetUpBC - Setup the subsections describing boundary conditions.
1251: Not Collective
1253: Input Parameter:
1254: . s - the `PetscSection`
1256: Level: advanced
1258: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSetUp()`, `PetscSectionCreate()`
1259: @*/
1260: PetscErrorCode PetscSectionSetUpBC(PetscSection s)
1261: {
1262: PetscFunctionBegin;
1264: if (s->bc) {
1265: const PetscInt last = (s->bc->pEnd - s->bc->pStart) - 1;
1267: PetscCall(PetscSectionSetUp(s->bc));
1268: if (last >= 0) PetscCall(PetscMalloc1(s->bc->atlasOff[last] + s->bc->atlasDof[last], &s->bcIndices));
1269: else s->bcIndices = NULL;
1270: }
1271: PetscFunctionReturn(PETSC_SUCCESS);
1272: }
1274: /*@
1275: PetscSectionSetUp - Calculate offsets based upon the number of degrees of freedom for each point in preparation for use of the `PetscSection`
1277: Not Collective
1279: Input Parameter:
1280: . s - the `PetscSection`
1282: Level: intermediate
1284: Notes:
1285: If used, `PetscSectionSetPermutation()` must be called before this routine.
1287: `PetscSectionSetPointMajor()`, cannot be called after this routine.
1289: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionSetPermutation()`
1290: @*/
1291: PetscErrorCode PetscSectionSetUp(PetscSection s)
1292: {
1293: PetscInt f;
1294: const PetscInt *pind = NULL;
1295: PetscCount offset = 0;
1297: PetscFunctionBegin;
1299: if (s->setup) PetscFunctionReturn(PETSC_SUCCESS);
1300: s->setup = PETSC_TRUE;
1301: /* Set offsets and field offsets for all points */
1302: /* Assume that all fields have the same chart */
1303: PetscCheck(s->includesConstraints, PETSC_COMM_SELF, PETSC_ERR_SUP, "PetscSectionSetUp is currently unsupported for includesConstraints = PETSC_TRUE");
1304: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1305: if (s->pointMajor) {
1306: PetscCount foff;
1307: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1308: const PetscInt q = pind ? pind[p] : p;
1310: /* Set point offset */
1311: PetscCall(PetscIntCast(offset, &s->atlasOff[q]));
1312: offset += s->atlasDof[q];
1313: /* Set field offset */
1314: for (f = 0, foff = s->atlasOff[q]; f < s->numFields; ++f) {
1315: PetscSection sf = s->field[f];
1317: PetscCall(PetscIntCast(foff, &sf->atlasOff[q]));
1318: foff += sf->atlasDof[q];
1319: }
1320: }
1321: } else {
1322: /* Set field offsets for all points */
1323: for (f = 0; f < s->numFields; ++f) {
1324: PetscSection sf = s->field[f];
1326: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1327: const PetscInt q = pind ? pind[p] : p;
1329: PetscCall(PetscIntCast(offset, &sf->atlasOff[q]));
1330: offset += sf->atlasDof[q];
1331: }
1332: }
1333: /* Disable point offsets since these are unused */
1334: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) s->atlasOff[p] = -1;
1335: }
1336: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1337: /* Setup BC sections */
1338: PetscCall(PetscSectionSetUpBC(s));
1339: for (f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetUpBC(s->field[f]));
1340: PetscFunctionReturn(PETSC_SUCCESS);
1341: }
1343: /*@
1344: PetscSectionGetMaxDof - Return the maximum number of degrees of freedom on any point in the `PetscSection`
1346: Not Collective
1348: Input Parameter:
1349: . s - the `PetscSection`
1351: Output Parameter:
1352: . maxDof - the maximum dof
1354: Level: intermediate
1356: Notes:
1357: The returned number is up-to-date without need for `PetscSectionSetUp()`.
1359: 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
1360: the maximum over all points of the value returned by `PetscSectionGetDof()` on this MPI process
1362: Developer Notes:
1363: The returned number is calculated lazily and stashed.
1365: A call to `PetscSectionInvalidateMaxDof_Internal()` invalidates the stashed value.
1367: `PetscSectionInvalidateMaxDof_Internal()` is called in `PetscSectionSetDof()`, `PetscSectionAddDof()` and `PetscSectionReset()`
1369: It should also be called every time `atlasDof` is modified directly.
1371: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
1372: @*/
1373: PetscErrorCode PetscSectionGetMaxDof(PetscSection s, PetscInt *maxDof)
1374: {
1375: PetscInt p;
1377: PetscFunctionBegin;
1379: PetscAssertPointer(maxDof, 2);
1380: if (s->maxDof == PETSC_INT_MIN) {
1381: s->maxDof = 0;
1382: for (p = 0; p < s->pEnd - s->pStart; ++p) s->maxDof = PetscMax(s->maxDof, s->atlasDof[p]);
1383: }
1384: *maxDof = s->maxDof;
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: /*@
1389: PetscSectionGetStorageSize - Return the size of an array or local `Vec` capable of holding all the degrees of freedom defined in a `PetscSection`
1391: Not Collective
1393: Input Parameter:
1394: . s - the `PetscSection`
1396: Output Parameter:
1397: . size - the size of an array which can hold all the dofs
1399: Level: intermediate
1401: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionGetConstrainedStorageSize()`, `PetscSectionCreate()`
1402: @*/
1403: PetscErrorCode PetscSectionGetStorageSize(PetscSection s, PetscInt *size)
1404: {
1405: PetscInt64 n = 0;
1407: PetscFunctionBegin;
1409: PetscAssertPointer(size, 2);
1410: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) n += s->atlasDof[p] > 0 ? s->atlasDof[p] : 0;
1411: PetscCall(PetscIntCast(n, size));
1412: PetscFunctionReturn(PETSC_SUCCESS);
1413: }
1415: /*@
1416: PetscSectionGetConstrainedStorageSize - Return the size of an array or local `Vec` capable of holding all unconstrained degrees of freedom in a `PetscSection`
1418: Not Collective
1420: Input Parameter:
1421: . s - the `PetscSection`
1423: Output Parameter:
1424: . size - the size of an array which can hold all unconstrained dofs
1426: Level: intermediate
1428: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetStorageSize()`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1429: @*/
1430: PetscErrorCode PetscSectionGetConstrainedStorageSize(PetscSection s, PetscInt *size)
1431: {
1432: PetscInt64 n = 0;
1434: PetscFunctionBegin;
1436: PetscAssertPointer(size, 2);
1437: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
1438: const PetscInt cdof = s->bc ? s->bc->atlasDof[p] : 0;
1439: n += s->atlasDof[p] > 0 ? s->atlasDof[p] - cdof : 0;
1440: }
1441: PetscCall(PetscIntCast(n, size));
1442: PetscFunctionReturn(PETSC_SUCCESS);
1443: }
1445: /*@
1446: PetscSectionCreateGlobalSection - Create a parallel section describing the global layout using
1447: a local (sequential) `PetscSection` on each MPI process and a `PetscSF` describing the section point overlap.
1449: Input Parameters:
1450: + s - The `PetscSection` for the local field layout
1451: . sf - The `PetscSF` describing parallel layout of the section points (leaves are unowned local points)
1452: . usePermutation - By default this is `PETSC_TRUE`, meaning any permutation of the local section is transferred to the global section
1453: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global field vector will not possess constrained dofs
1454: - localOffsets - If `PETSC_TRUE`, use local rather than global offsets for the points
1456: Output Parameter:
1457: . gsection - The `PetscSection` for the global field layout
1459: Level: intermediate
1461: Notes:
1462: On each MPI process `gsection` inherits the chart of the `s` on that process.
1464: 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`.
1465: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1467: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCreateGlobalSectionCensored()`
1468: @*/
1469: PetscErrorCode PetscSectionCreateGlobalSection(PetscSection s, PetscSF sf, PetscBool usePermutation, PetscBool includeConstraints, PetscBool localOffsets, PetscSection *gsection)
1470: {
1471: PetscSection gs;
1472: const PetscInt *pind = NULL;
1473: PetscInt *recv = NULL, *neg = NULL;
1474: PetscInt pStart, pEnd, p, dof, cdof, off, globalOff = 0, nroots, nlocal, maxleaf;
1475: PetscInt numFields, f, numComponents;
1476: PetscCount foff;
1478: PetscFunctionBegin;
1484: PetscAssertPointer(gsection, 6);
1485: PetscCheck(s->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
1486: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &gs));
1487: PetscCall(PetscSectionGetNumFields(s, &numFields));
1488: if (numFields > 0) PetscCall(PetscSectionSetNumFields(gs, numFields));
1489: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1490: PetscCall(PetscSectionSetChart(gs, pStart, pEnd));
1491: gs->includesConstraints = includeConstraints;
1492: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1493: nlocal = nroots; /* The local/leaf space matches global/root space */
1494: /* Must allocate for all points visible to SF, which may be more than this section */
1495: if (nroots >= 0) { /* nroots < 0 means that the graph has not been set, only happens in serial */
1496: PetscCall(PetscSFGetLeafRange(sf, NULL, &maxleaf));
1497: PetscCheck(nroots >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "SF roots %" PetscInt_FMT " < pEnd %" PetscInt_FMT, nroots, pEnd);
1498: PetscCheck(maxleaf < nroots, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Max local leaf %" PetscInt_FMT " >= nroots %" PetscInt_FMT, maxleaf, nroots);
1499: PetscCall(PetscMalloc2(nroots, &neg, nlocal, &recv));
1500: PetscCall(PetscArrayzero(neg, nroots));
1501: }
1502: /* Mark all local points with negative dof */
1503: for (p = pStart; p < pEnd; ++p) {
1504: PetscCall(PetscSectionGetDof(s, p, &dof));
1505: PetscCall(PetscSectionSetDof(gs, p, dof));
1506: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1507: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(gs, p, cdof));
1508: if (neg) neg[p] = -(dof + 1);
1509: }
1510: PetscCall(PetscSectionSetUpBC(gs));
1511: 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]));
1512: if (nroots >= 0) {
1513: PetscCall(PetscArrayzero(recv, nlocal));
1514: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1515: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1516: for (p = pStart; p < pEnd; ++p) {
1517: if (recv[p] < 0) {
1518: gs->atlasDof[p - pStart] = recv[p];
1519: PetscCall(PetscSectionGetDof(s, p, &dof));
1520: 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);
1521: }
1522: }
1523: }
1524: /* Calculate new sizes, get process offset, and calculate point offsets */
1525: if (usePermutation && s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1526: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1527: const PetscInt q = pind ? pind[p] : p;
1529: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1530: gs->atlasOff[q] = off;
1531: off += gs->atlasDof[q] > 0 ? gs->atlasDof[q] - cdof : 0;
1532: }
1533: if (!localOffsets) {
1534: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1535: globalOff -= off;
1536: }
1537: for (p = pStart, off = 0; p < pEnd; ++p) {
1538: gs->atlasOff[p - pStart] += globalOff;
1539: if (neg) neg[p] = -(gs->atlasOff[p - pStart] + 1);
1540: }
1541: if (usePermutation && s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1542: /* Put in negative offsets for ghost points */
1543: if (nroots >= 0) {
1544: PetscCall(PetscArrayzero(recv, nlocal));
1545: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1546: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, recv, MPI_REPLACE));
1547: for (p = pStart; p < pEnd; ++p) {
1548: if (recv[p] < 0) gs->atlasOff[p - pStart] = recv[p];
1549: }
1550: }
1551: PetscCall(PetscFree2(neg, recv));
1552: /* Set field dofs/offsets/constraints */
1553: for (f = 0; f < numFields; ++f) {
1554: const char *name;
1556: gs->field[f]->includesConstraints = includeConstraints;
1557: PetscCall(PetscSectionGetFieldComponents(s, f, &numComponents));
1558: PetscCall(PetscSectionSetFieldComponents(gs, f, numComponents));
1559: PetscCall(PetscSectionGetFieldName(s, f, &name));
1560: PetscCall(PetscSectionSetFieldName(gs, f, name));
1561: }
1562: for (p = pStart; p < pEnd; ++p) {
1563: PetscCall(PetscSectionGetOffset(gs, p, &off));
1564: for (f = 0, foff = off; f < numFields; ++f) {
1565: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
1566: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetFieldConstraintDof(gs, p, f, cdof));
1567: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
1568: PetscCall(PetscSectionSetFieldDof(gs, p, f, off < 0 ? -(dof + 1) : dof));
1569: PetscCall(PetscSectionSetFieldOffset(gs, p, f, (PetscInt)foff));
1570: PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof));
1571: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1572: PetscCheck(foff < PETSC_INT_MAX, PETSC_COMM_SELF, PETSC_ERR_INT_OVERFLOW, "Offsets too large for 32 bit indices");
1573: }
1574: }
1575: for (f = 0; f < numFields; ++f) {
1576: PetscSection gfs = gs->field[f];
1578: PetscCall(PetscSectionSetUpBC(gfs));
1579: 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]));
1580: }
1581: gs->setup = PETSC_TRUE;
1582: PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view"));
1583: *gsection = gs;
1584: PetscFunctionReturn(PETSC_SUCCESS);
1585: }
1587: /*@
1588: PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the globallayout using
1589: a local (sequential) `PetscSection` on each MPI process and an `PetscSF` describing the section point overlap.
1591: Input Parameters:
1592: + s - The `PetscSection` for the local field layout
1593: . sf - The `PetscSF` describing parallel layout of the section points
1594: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global vector will not possess constrained dofs
1595: . numExcludes - The number of exclusion ranges, this must have the same value on all MPI processes
1596: - 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
1598: Output Parameter:
1599: . gsection - The `PetscSection` for the global field layout
1601: Level: advanced
1603: Notes:
1604: On each MPI process `gsection` inherits the chart of the `s` on that process.
1606: 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`.
1607: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1609: This routine augments `PetscSectionCreateGlobalSection()` by allowing one to exclude certain ranges in the chart of the `PetscSection`
1611: Developer Notes:
1612: This is a terrible function name
1614: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
1615: @*/
1616: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1617: {
1618: const PetscInt *pind = NULL;
1619: PetscInt *neg = NULL, *tmpOff = NULL;
1620: PetscInt pStart, pEnd, p, e, dof, cdof, globalOff = 0, nroots;
1621: PetscCount off;
1623: PetscFunctionBegin;
1626: PetscAssertPointer(gsection, 6);
1627: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
1628: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1629: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
1630: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1631: if (nroots >= 0) {
1632: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
1633: PetscCall(PetscCalloc1(nroots, &neg));
1634: if (nroots > pEnd - pStart) {
1635: PetscCall(PetscCalloc1(nroots, &tmpOff));
1636: } else {
1637: tmpOff = &(*gsection)->atlasDof[-pStart];
1638: }
1639: }
1640: /* Mark ghost points with negative dof */
1641: for (p = pStart; p < pEnd; ++p) {
1642: for (e = 0; e < numExcludes; ++e) {
1643: if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1644: PetscCall(PetscSectionSetDof(*gsection, p, 0));
1645: break;
1646: }
1647: }
1648: if (e < numExcludes) continue;
1649: PetscCall(PetscSectionGetDof(s, p, &dof));
1650: PetscCall(PetscSectionSetDof(*gsection, p, dof));
1651: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1652: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
1653: if (neg) neg[p] = -(dof + 1);
1654: }
1655: PetscCall(PetscSectionSetUpBC(*gsection));
1656: if (nroots >= 0) {
1657: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1658: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1659: if (nroots > pEnd - pStart) {
1660: for (p = pStart; p < pEnd; ++p) {
1661: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1662: }
1663: }
1664: }
1665: /* Calculate new sizes, get process offset, and calculate point offsets */
1666: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1667: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1668: const PetscInt q = pind ? pind[p] : p;
1670: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1671: (*gsection)->atlasOff[q] = (PetscInt)off;
1672: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1673: PetscCheck(off < PETSC_INT_MAX, PETSC_COMM_SELF, PETSC_ERR_INT_OVERFLOW, "Offsets too large for 32 bit indices");
1674: }
1675: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1676: globalOff -= off;
1677: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1678: (*gsection)->atlasOff[p] += globalOff;
1679: if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1680: }
1681: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1682: /* Put in negative offsets for ghost points */
1683: if (nroots >= 0) {
1684: if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1685: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1686: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1687: if (nroots > pEnd - pStart) {
1688: for (p = pStart; p < pEnd; ++p) {
1689: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1690: }
1691: }
1692: }
1693: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
1694: PetscCall(PetscFree(neg));
1695: PetscFunctionReturn(PETSC_SUCCESS);
1696: }
1698: /*@
1699: PetscSectionGetPointLayout - Get a `PetscLayout` for the points with nonzero dof counts of the unnamed default field within this `PetscSection`s local chart
1701: Collective
1703: Input Parameters:
1704: + comm - The `MPI_Comm`
1705: - s - The `PetscSection`
1707: Output Parameter:
1708: . layout - The point layout for the data that defines the section
1710: Level: advanced
1712: Notes:
1713: `PetscSectionGetValueLayout()` provides similar information but counting the total number of degrees of freedom on the MPI process (excluding constrained
1714: degrees of freedom).
1716: This count includes constrained degrees of freedom
1718: This is usually called on the default global section.
1720: Example:
1721: .vb
1722: The chart is [2,5), point 2 has 2 dof, point 3 has 0 dof, point 4 has 1 dof
1723: The local size of the `PetscLayout` is 2 since 2 points have a non-zero number of dof
1724: .ve
1726: Developer Notes:
1727: I find the names of these two functions extremely non-informative
1729: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1730: @*/
1731: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1732: {
1733: PetscInt pStart, pEnd, p, localSize = 0;
1735: PetscFunctionBegin;
1736: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1737: for (p = pStart; p < pEnd; ++p) {
1738: PetscInt dof;
1740: PetscCall(PetscSectionGetDof(s, p, &dof));
1741: if (dof >= 0) ++localSize;
1742: }
1743: PetscCall(PetscLayoutCreate(comm, layout));
1744: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1745: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1746: PetscCall(PetscLayoutSetUp(*layout));
1747: PetscFunctionReturn(PETSC_SUCCESS);
1748: }
1750: /*@
1751: PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs of a `PetscSection`
1753: Collective
1755: Input Parameters:
1756: + comm - The `MPI_Comm`
1757: - s - The `PetscSection`
1759: Output Parameter:
1760: . layout - The dof layout for the section
1762: Level: advanced
1764: Notes:
1765: `PetscSectionGetPointLayout()` provides similar information but only counting the number of points with nonzero degrees of freedom and
1766: including the constrained degrees of freedom
1768: This is usually called for the default global section.
1770: Example:
1771: .vb
1772: The chart is [2,5), point 2 has 4 dof (2 constrained), point 3 has 0 dof, point 4 has 1 dof (not constrained)
1773: The local size of the `PetscLayout` is 3 since there are 3 unconstrained degrees of freedom on this MPI process
1774: .ve
1776: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1777: @*/
1778: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1779: {
1780: PetscInt pStart, pEnd, p, localSize = 0;
1782: PetscFunctionBegin;
1784: PetscAssertPointer(layout, 3);
1785: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1786: for (p = pStart; p < pEnd; ++p) {
1787: PetscInt dof, cdof;
1789: PetscCall(PetscSectionGetDof(s, p, &dof));
1790: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1791: if (dof - cdof > 0) localSize += dof - cdof;
1792: }
1793: PetscCall(PetscLayoutCreate(comm, layout));
1794: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1795: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1796: PetscCall(PetscLayoutSetUp(*layout));
1797: PetscFunctionReturn(PETSC_SUCCESS);
1798: }
1800: /*@
1801: PetscSectionGetOffset - Return the offset into an array or `Vec` for the dof associated with the given point.
1803: Not Collective
1805: Input Parameters:
1806: + s - the `PetscSection`
1807: - point - the point
1809: Output Parameter:
1810: . offset - the offset
1812: Level: intermediate
1814: Notes:
1815: In a global section, `offset` will be negative for points not owned by this process.
1817: This is for the unnamed default field in the `PetscSection` not the named fields
1819: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1821: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetPointMajor()`
1822: @*/
1823: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1824: {
1825: PetscFunctionBegin;
1827: PetscAssertPointer(offset, 3);
1828: 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);
1829: *offset = s->atlasOff[point - s->pStart];
1830: PetscFunctionReturn(PETSC_SUCCESS);
1831: }
1833: /*@
1834: PetscSectionSetOffset - Set the offset into an array or `Vec` for the dof associated with the given point.
1836: Not Collective
1838: Input Parameters:
1839: + s - the `PetscSection`
1840: . point - the point
1841: - offset - the offset, these values may be negative indicating the values are off process
1843: Level: developer
1845: Note:
1846: The user usually does not call this function, but uses `PetscSectionSetUp()`
1848: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1849: @*/
1850: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1851: {
1852: PetscFunctionBegin;
1854: 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);
1855: s->atlasOff[point - s->pStart] = offset;
1856: PetscFunctionReturn(PETSC_SUCCESS);
1857: }
1859: /*@
1860: PetscSectionGetFieldOffset - Return the offset into an array or `Vec` for the field dof associated with the given point.
1862: Not Collective
1864: Input Parameters:
1865: + s - the `PetscSection`
1866: . point - the point
1867: - field - the field
1869: Output Parameter:
1870: . offset - the offset
1872: Level: intermediate
1874: Notes:
1875: In a global section, `offset` will be negative for points not owned by this process.
1877: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1879: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldPointOffset()`
1880: @*/
1881: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1882: {
1883: PetscFunctionBegin;
1885: PetscAssertPointer(offset, 4);
1886: PetscSectionCheckValidField(field, s->numFields);
1887: PetscCall(PetscSectionGetOffset(s->field[field], point, offset));
1888: PetscFunctionReturn(PETSC_SUCCESS);
1889: }
1891: /*@
1892: PetscSectionSetFieldOffset - Set the offset into an array or `Vec` for the dof associated with the given field at a point.
1894: Not Collective
1896: Input Parameters:
1897: + s - the `PetscSection`
1898: . point - the point
1899: . field - the field
1900: - offset - the offset, these values may be negative indicating the values are off process
1902: Level: developer
1904: Note:
1905: The user usually does not call this function, but uses `PetscSectionSetUp()`
1907: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1908: @*/
1909: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1910: {
1911: PetscFunctionBegin;
1913: PetscSectionCheckValidField(field, s->numFields);
1914: PetscCall(PetscSectionSetOffset(s->field[field], point, offset));
1915: PetscFunctionReturn(PETSC_SUCCESS);
1916: }
1918: /*@
1919: PetscSectionGetFieldPointOffset - Return the offset for the first field dof associated with the given point relative to the offset for that point for the
1920: unnamed default field's first dof
1922: Not Collective
1924: Input Parameters:
1925: + s - the `PetscSection`
1926: . point - the point
1927: - field - the field
1929: Output Parameter:
1930: . offset - the offset
1932: Level: advanced
1934: Note:
1935: This ignores constraints
1937: Example:
1938: .vb
1939: if PetscSectionSetPointMajor(s,PETSC_TRUE)
1940: The unnamed default field has 3 dof at `point`
1941: Field 0 has 2 dof at `point`
1942: Then PetscSectionGetFieldPointOffset(s,point,1,&offset) returns and offset of 5
1943: .ve
1945: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldOffset()`
1946: @*/
1947: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1948: {
1949: PetscInt off, foff;
1951: PetscFunctionBegin;
1953: PetscAssertPointer(offset, 4);
1954: PetscSectionCheckValidField(field, s->numFields);
1955: PetscCall(PetscSectionGetOffset(s, point, &off));
1956: PetscCall(PetscSectionGetOffset(s->field[field], point, &foff));
1957: *offset = foff - off;
1958: PetscFunctionReturn(PETSC_SUCCESS);
1959: }
1961: /*@
1962: PetscSectionGetOffsetRange - Return the full range of offsets [`start`, `end`) for a `PetscSection`
1964: Not Collective
1966: Input Parameter:
1967: . s - the `PetscSection`
1969: Output Parameters:
1970: + start - the minimum offset
1971: - end - one more than the maximum offset
1973: Level: intermediate
1975: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1976: @*/
1977: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1978: {
1979: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1981: PetscFunctionBegin;
1983: if (s->atlasOff) {
1984: os = s->atlasOff[0];
1985: oe = s->atlasOff[0];
1986: }
1987: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1988: for (p = 0; p < pEnd - pStart; ++p) {
1989: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1991: if (off >= 0) {
1992: os = PetscMin(os, off);
1993: oe = PetscMax(oe, off + dof);
1994: }
1995: }
1996: if (start) *start = os;
1997: if (end) *end = oe;
1998: PetscFunctionReturn(PETSC_SUCCESS);
1999: }
2001: /*@
2002: PetscSectionCreateSubsection - Create a new, smaller `PetscSection` composed of only selected fields
2004: Collective
2006: Input Parameters:
2007: + s - the `PetscSection`
2008: . len - the number of subfields
2009: - fields - the subfield numbers
2011: Output Parameter:
2012: . subs - the subsection
2014: Level: advanced
2016: Notes:
2017: The chart of `subs` is the same as the chart of `s`
2019: This will error if a fieldnumber is out of range
2021: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2022: @*/
2023: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
2024: {
2025: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
2027: PetscFunctionBegin;
2028: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2030: PetscAssertPointer(fields, 3);
2031: PetscAssertPointer(subs, 4);
2032: PetscCall(PetscSectionGetNumFields(s, &nF));
2033: 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);
2034: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2035: PetscCall(PetscSectionSetNumFields(*subs, len));
2036: for (f = 0; f < len; ++f) {
2037: const char *name = NULL;
2038: PetscInt numComp = 0;
2039: PetscSectionSym sym;
2041: PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
2042: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2043: PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
2044: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2045: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
2046: PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
2047: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2048: }
2049: PetscCall(PetscSectionGetFieldSym(s, fields[f], &sym));
2050: PetscCall(PetscSectionSetFieldSym(*subs, f, sym));
2051: }
2052: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2053: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2054: for (p = pStart; p < pEnd; ++p) {
2055: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
2057: for (f = 0; f < len; ++f) {
2058: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2059: PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
2060: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2061: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
2062: dof += fdof;
2063: cdof += cfdof;
2064: }
2065: PetscCall(PetscSectionSetDof(*subs, p, dof));
2066: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
2067: maxCdof = PetscMax(cdof, maxCdof);
2068: }
2069: PetscCall(PetscSectionSetUp(*subs));
2070: if (maxCdof) {
2071: PetscInt *indices;
2073: PetscCall(PetscMalloc1(maxCdof, &indices));
2074: for (p = pStart; p < pEnd; ++p) {
2075: PetscInt cdof;
2077: PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
2078: if (cdof) {
2079: const PetscInt *oldIndices = NULL;
2080: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
2082: for (f = 0; f < len; ++f) {
2083: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2084: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2085: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
2086: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
2087: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
2088: numConst += cfdof;
2089: fOff += fdof;
2090: }
2091: PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
2092: }
2093: }
2094: PetscCall(PetscFree(indices));
2095: }
2096: PetscFunctionReturn(PETSC_SUCCESS);
2097: }
2099: /*@
2100: PetscSectionCreateComponentSubsection - Create a new, smaller `PetscSection` composed of only selected components
2102: Collective
2104: Input Parameters:
2105: + s - the `PetscSection`
2106: . len - the number of components
2107: - comps - the component numbers
2109: Output Parameter:
2110: . subs - the subsection
2112: Level: advanced
2114: Notes:
2115: The chart of `subs` is the same as the chart of `s`
2117: This will error if the section has more than one field, or if a component number is out of range
2119: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2120: @*/
2121: PetscErrorCode PetscSectionCreateComponentSubsection(PetscSection s, PetscInt len, const PetscInt comps[], PetscSection *subs)
2122: {
2123: PetscSectionSym sym;
2124: const char *name = NULL;
2125: PetscInt Nf, pStart, pEnd;
2127: PetscFunctionBegin;
2128: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2130: PetscAssertPointer(comps, 3);
2131: PetscAssertPointer(subs, 4);
2132: PetscCall(PetscSectionGetNumFields(s, &Nf));
2133: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "This method can only handle one field, not %" PetscInt_FMT, Nf);
2134: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2135: PetscCall(PetscSectionSetNumFields(*subs, 1));
2136: PetscCall(PetscSectionGetFieldName(s, 0, &name));
2137: PetscCall(PetscSectionSetFieldName(*subs, 0, name));
2138: PetscCall(PetscSectionSetFieldComponents(*subs, 0, len));
2139: PetscCall(PetscSectionGetFieldSym(s, 0, &sym));
2140: PetscCall(PetscSectionSetFieldSym(*subs, 0, sym));
2141: for (PetscInt c = 0; c < len; ++c) {
2142: PetscCall(PetscSectionGetComponentName(s, 0, comps[c], &name));
2143: PetscCall(PetscSectionSetComponentName(*subs, 0, c, name));
2144: }
2145: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2146: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2147: for (PetscInt p = pStart; p < pEnd; ++p) {
2148: PetscInt dof, cdof, cfdof;
2150: PetscCall(PetscSectionGetDof(s, p, &dof));
2151: if (!dof) continue;
2152: PetscCall(PetscSectionGetFieldConstraintDof(s, p, 0, &cfdof));
2153: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2154: PetscCheck(!cdof && !cfdof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Component selection does not work with constraints");
2155: PetscCall(PetscSectionSetFieldDof(*subs, p, 0, len));
2156: PetscCall(PetscSectionSetDof(*subs, p, len));
2157: }
2158: PetscCall(PetscSectionSetUp(*subs));
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: /*@
2163: PetscSectionCreateSupersection - Create a new, larger section composed of multiple `PetscSection`s
2165: Collective
2167: Input Parameters:
2168: + s - the input sections
2169: - len - the number of input sections
2171: Output Parameter:
2172: . supers - the supersection
2174: Level: advanced
2176: Notes:
2177: The section offsets now refer to a new, larger vector.
2179: Developer Notes:
2180: Needs to explain how the sections are composed
2182: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2183: @*/
2184: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2185: {
2186: PetscInt Nf = 0, f, pStart = PETSC_INT_MAX, pEnd = 0, p, maxCdof = 0, i;
2188: PetscFunctionBegin;
2189: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2190: for (i = 0; i < len; ++i) {
2191: PetscInt nf, pStarti, pEndi;
2193: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2194: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2195: pStart = PetscMin(pStart, pStarti);
2196: pEnd = PetscMax(pEnd, pEndi);
2197: Nf += nf;
2198: }
2199: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2200: PetscCall(PetscSectionSetNumFields(*supers, Nf));
2201: for (i = 0, f = 0; i < len; ++i) {
2202: PetscInt nf, fi, ci;
2204: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2205: for (fi = 0; fi < nf; ++fi, ++f) {
2206: const char *name = NULL;
2207: PetscInt numComp = 0;
2209: PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2210: PetscCall(PetscSectionSetFieldName(*supers, f, name));
2211: PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2212: PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2213: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2214: PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2215: PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2216: }
2217: }
2218: }
2219: PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2220: for (p = pStart; p < pEnd; ++p) {
2221: PetscInt dof = 0, cdof = 0;
2223: for (i = 0, f = 0; i < len; ++i) {
2224: PetscInt nf, fi, pStarti, pEndi;
2225: PetscInt fdof = 0, cfdof = 0;
2227: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2228: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2229: if ((p < pStarti) || (p >= pEndi)) continue;
2230: for (fi = 0; fi < nf; ++fi, ++f) {
2231: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2232: PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2233: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2234: if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2235: dof += fdof;
2236: cdof += cfdof;
2237: }
2238: }
2239: PetscCall(PetscSectionSetDof(*supers, p, dof));
2240: if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2241: maxCdof = PetscMax(cdof, maxCdof);
2242: }
2243: PetscCall(PetscSectionSetUp(*supers));
2244: if (maxCdof) {
2245: PetscInt *indices;
2247: PetscCall(PetscMalloc1(maxCdof, &indices));
2248: for (p = pStart; p < pEnd; ++p) {
2249: PetscInt cdof;
2251: PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2252: if (cdof) {
2253: PetscInt dof, numConst = 0, fOff = 0;
2255: for (i = 0, f = 0; i < len; ++i) {
2256: const PetscInt *oldIndices = NULL;
2257: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
2259: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2260: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2261: if ((p < pStarti) || (p >= pEndi)) continue;
2262: for (fi = 0; fi < nf; ++fi, ++f) {
2263: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2264: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2265: PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2266: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2267: PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2268: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2269: numConst += cfdof;
2270: }
2271: PetscCall(PetscSectionGetDof(s[i], p, &dof));
2272: fOff += dof;
2273: }
2274: PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2275: }
2276: }
2277: PetscCall(PetscFree(indices));
2278: }
2279: PetscFunctionReturn(PETSC_SUCCESS);
2280: }
2282: static PetscErrorCode PetscSectionCreateSubplexSection_Private(PetscSection s, IS subpointIS, PetscBool renumberPoints, PetscSection *subs)
2283: {
2284: const PetscInt *points = NULL, *indices = NULL;
2285: PetscInt *spoints = NULL, *order = NULL;
2286: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
2288: PetscFunctionBegin;
2291: PetscAssertPointer(subs, 4);
2292: PetscCall(PetscSectionGetNumFields(s, &numFields));
2293: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2294: if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2295: for (f = 0; f < numFields; ++f) {
2296: const char *name = NULL;
2297: PetscInt numComp = 0;
2299: PetscCall(PetscSectionGetFieldName(s, f, &name));
2300: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2301: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2302: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2303: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2304: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2305: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2306: }
2307: }
2308: /* For right now, we do not try to squeeze the subchart */
2309: if (subpointIS) {
2310: PetscCall(ISGetLocalSize(subpointIS, &numSubpoints));
2311: PetscCall(ISGetIndices(subpointIS, &points));
2312: }
2313: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2314: if (renumberPoints) {
2315: PetscBool sorted;
2317: spStart = 0;
2318: spEnd = numSubpoints;
2319: PetscCall(ISSorted(subpointIS, &sorted));
2320: if (!sorted) {
2321: PetscCall(PetscMalloc2(numSubpoints, &spoints, numSubpoints, &order));
2322: PetscCall(PetscArraycpy(spoints, points, numSubpoints));
2323: for (PetscInt i = 0; i < numSubpoints; ++i) order[i] = i;
2324: PetscCall(PetscSortIntWithArray(numSubpoints, spoints, order));
2325: }
2326: } else {
2327: PetscCall(ISGetMinMax(subpointIS, &spStart, &spEnd));
2328: ++spEnd;
2329: }
2330: PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2331: for (p = pStart; p < pEnd; ++p) {
2332: PetscInt dof, cdof, fdof = 0, cfdof = 0;
2334: PetscCall(PetscFindInt(p, numSubpoints, spoints ? spoints : points, &subp));
2335: if (subp < 0) continue;
2336: if (!renumberPoints) subp = p;
2337: else subp = order ? order[subp] : subp;
2338: for (f = 0; f < numFields; ++f) {
2339: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2340: PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2341: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2342: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2343: }
2344: PetscCall(PetscSectionGetDof(s, p, &dof));
2345: PetscCall(PetscSectionSetDof(*subs, subp, dof));
2346: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2347: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2348: }
2349: PetscCall(PetscSectionSetUp(*subs));
2350: /* Change offsets to original offsets */
2351: for (p = pStart; p < pEnd; ++p) {
2352: PetscInt off, foff = 0;
2354: PetscCall(PetscFindInt(p, numSubpoints, spoints ? spoints : points, &subp));
2355: if (subp < 0) continue;
2356: if (!renumberPoints) subp = p;
2357: else subp = order ? order[subp] : subp;
2358: for (f = 0; f < numFields; ++f) {
2359: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2360: PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2361: }
2362: PetscCall(PetscSectionGetOffset(s, p, &off));
2363: PetscCall(PetscSectionSetOffset(*subs, subp, off));
2364: }
2365: /* Copy constraint indices */
2366: for (subp = spStart; subp < spEnd; ++subp) {
2367: PetscInt cdof;
2369: PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2370: if (cdof) {
2371: for (f = 0; f < numFields; ++f) {
2372: PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2373: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2374: }
2375: PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2376: PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2377: }
2378: }
2379: if (subpointIS) PetscCall(ISRestoreIndices(subpointIS, &points));
2380: PetscCall(PetscFree2(spoints, order));
2381: PetscFunctionReturn(PETSC_SUCCESS);
2382: }
2384: /*@
2385: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2387: Collective
2389: Input Parameters:
2390: + s - the `PetscSection`
2391: - subpointIS - a sorted list of points in the original mesh which are in the submesh
2393: Output Parameter:
2394: . subs - the subsection
2396: Level: advanced
2398: Notes:
2399: 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))`
2401: Compare this with `PetscSectionCreateSubdomainSection()` that does not map the points numbers to start at zero but leaves them as before
2403: Developer Notes:
2404: 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`
2406: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2407: @*/
2408: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointIS, PetscSection *subs)
2409: {
2410: PetscFunctionBegin;
2411: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointIS, PETSC_TRUE, subs));
2412: PetscFunctionReturn(PETSC_SUCCESS);
2413: }
2415: /*@
2416: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2418: Collective
2420: Input Parameters:
2421: + s - the `PetscSection`
2422: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2424: Output Parameter:
2425: . subs - the subsection
2427: Level: advanced
2429: Notes:
2430: 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`
2431: is `[min(subpointMap),max(subpointMap)+1)`
2433: Compare this with `PetscSectionCreateSubmeshSection()` that maps the point numbers to start at zero
2435: Developer Notes:
2436: 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`
2438: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2439: @*/
2440: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2441: {
2442: PetscFunctionBegin;
2443: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_FALSE, subs));
2444: PetscFunctionReturn(PETSC_SUCCESS);
2445: }
2447: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2448: {
2449: PetscInt p;
2450: PetscMPIInt rank;
2452: PetscFunctionBegin;
2453: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2454: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2455: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2456: for (p = 0; p < s->pEnd - s->pStart; ++p) {
2457: if (s->bc && s->bc->atlasDof[p] > 0) {
2458: PetscInt b;
2459: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2460: if (s->bcIndices) {
2461: for (b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2462: }
2463: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2464: } else {
2465: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2466: }
2467: }
2468: PetscCall(PetscViewerFlush(viewer));
2469: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2470: if (s->sym) {
2471: PetscCall(PetscViewerASCIIPushTab(viewer));
2472: PetscCall(PetscSectionSymView(s->sym, viewer));
2473: PetscCall(PetscViewerASCIIPopTab(viewer));
2474: }
2475: PetscFunctionReturn(PETSC_SUCCESS);
2476: }
2478: /*@
2479: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2481: Collective
2483: Input Parameters:
2484: + A - the `PetscSection` object to view
2485: . obj - Optional object that provides the options prefix used for the options
2486: - name - command line option
2488: Level: intermediate
2490: Note:
2491: See `PetscObjectViewFromOptions()` for available values of `PetscViewer` and `PetscViewerFormat`
2493: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2494: @*/
2495: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2496: {
2497: PetscFunctionBegin;
2499: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2500: PetscFunctionReturn(PETSC_SUCCESS);
2501: }
2503: /*@
2504: PetscSectionView - Views a `PetscSection`
2506: Collective
2508: Input Parameters:
2509: + s - the `PetscSection` object to view
2510: - viewer - the viewer
2512: Level: beginner
2514: Note:
2515: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2516: distribution independent data, such as dofs, offsets, constraint dofs,
2517: and constraint indices. Points that have negative dofs, for instance,
2518: are not saved as they represent points owned by other processes.
2519: Point numbering and rank assignment is currently not stored.
2520: The saved section can be loaded with `PetscSectionLoad()`.
2522: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2523: @*/
2524: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2525: {
2526: PetscBool isascii, ishdf5;
2527: PetscInt f;
2529: PetscFunctionBegin;
2531: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2533: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2534: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2535: if (isascii) {
2536: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2537: if (s->numFields) {
2538: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2539: for (f = 0; f < s->numFields; ++f) {
2540: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " \"%s\" with %" PetscInt_FMT " components\n", f, s->fieldNames[f], s->numFieldComponents[f]));
2541: PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2542: }
2543: } else {
2544: PetscCall(PetscSectionView_ASCII(s, viewer));
2545: }
2546: } else if (ishdf5) {
2547: #if PetscDefined(HAVE_HDF5)
2548: PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2549: #else
2550: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2551: #endif
2552: }
2553: PetscFunctionReturn(PETSC_SUCCESS);
2554: }
2556: /*@
2557: PetscSectionLoad - Loads a `PetscSection`
2559: Collective
2561: Input Parameters:
2562: + s - the `PetscSection` object to load
2563: - viewer - the viewer
2565: Level: beginner
2567: Note:
2568: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2569: a section saved with `PetscSectionView()`. The number of processes
2570: used here (N) does not need to be the same as that used when saving.
2571: After calling this function, the chart of s on rank i will be set
2572: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2573: saved section points.
2575: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2576: @*/
2577: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2578: {
2579: PetscBool ishdf5;
2581: PetscFunctionBegin;
2584: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2585: if (ishdf5) {
2586: #if PetscDefined(HAVE_HDF5)
2587: PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2588: PetscFunctionReturn(PETSC_SUCCESS);
2589: #else
2590: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2591: #endif
2592: } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2593: }
2595: /*@
2596: PetscSectionResetClosurePermutation - Remove any existing closure permutation
2598: Input Parameter:
2599: . section - The `PetscSection`
2601: Level: intermediate
2603: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionSetClosureIndex()`, `PetscSectionReset()`
2604: @*/
2605: PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2606: {
2607: PetscSectionClosurePermVal clVal;
2609: PetscFunctionBegin;
2610: if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2611: kh_foreach_value(section->clHash, clVal, {
2612: PetscCall(PetscFree(clVal.perm));
2613: PetscCall(PetscFree(clVal.invPerm));
2614: });
2615: kh_destroy(ClPerm, section->clHash);
2616: section->clHash = NULL;
2617: PetscFunctionReturn(PETSC_SUCCESS);
2618: }
2620: /*@
2621: PetscSectionReset - Frees all section data, the section is then as if `PetscSectionCreate()` had just been called.
2623: Not Collective
2625: Input Parameter:
2626: . s - the `PetscSection`
2628: Level: beginner
2630: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`
2631: @*/
2632: PetscErrorCode PetscSectionReset(PetscSection s)
2633: {
2634: PetscInt f, c;
2636: PetscFunctionBegin;
2638: for (f = 0; f < s->numFields; ++f) {
2639: PetscCall(PetscSectionDestroy(&s->field[f]));
2640: PetscCall(PetscFree(s->fieldNames[f]));
2641: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2642: PetscCall(PetscFree(s->compNames[f]));
2643: }
2644: PetscCall(PetscFree(s->numFieldComponents));
2645: PetscCall(PetscFree(s->fieldNames));
2646: PetscCall(PetscFree(s->compNames));
2647: PetscCall(PetscFree(s->field));
2648: PetscCall(PetscSectionDestroy(&s->bc));
2649: PetscCall(PetscFree(s->bcIndices));
2650: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2651: PetscCall(PetscSectionDestroy(&s->clSection));
2652: PetscCall(ISDestroy(&s->clPoints));
2653: PetscCall(ISDestroy(&s->perm));
2654: PetscCall(PetscBTDestroy(&s->blockStarts));
2655: PetscCall(PetscSectionResetClosurePermutation(s));
2656: PetscCall(PetscSectionSymDestroy(&s->sym));
2657: PetscCall(PetscSectionDestroy(&s->clSection));
2658: PetscCall(ISDestroy(&s->clPoints));
2659: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2660: s->pStart = -1;
2661: s->pEnd = -1;
2662: s->maxDof = 0;
2663: s->setup = PETSC_FALSE;
2664: s->numFields = 0;
2665: s->clObj = NULL;
2666: PetscFunctionReturn(PETSC_SUCCESS);
2667: }
2669: /*@
2670: PetscSectionDestroy - Frees a `PetscSection`
2672: Not Collective
2674: Input Parameter:
2675: . s - the `PetscSection`
2677: Level: beginner
2679: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2680: @*/
2681: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2682: {
2683: PetscFunctionBegin;
2684: if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2686: if (--((PetscObject)*s)->refct > 0) {
2687: *s = NULL;
2688: PetscFunctionReturn(PETSC_SUCCESS);
2689: }
2690: PetscCall(PetscSectionReset(*s));
2691: PetscCall(PetscHeaderDestroy(s));
2692: PetscFunctionReturn(PETSC_SUCCESS);
2693: }
2695: static PetscErrorCode VecIntGetValuesSection_Private(const PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2696: {
2697: const PetscInt p = point - s->pStart;
2699: PetscFunctionBegin;
2701: *values = &baseArray[s->atlasOff[p]];
2702: PetscFunctionReturn(PETSC_SUCCESS);
2703: }
2705: static PetscErrorCode VecIntSetValuesSection_Private(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2706: {
2707: PetscInt *array;
2708: const PetscInt p = point - s->pStart;
2709: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2710: PetscInt cDim = 0;
2712: PetscFunctionBegin;
2714: PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2715: array = &baseArray[s->atlasOff[p]];
2716: if (!cDim) {
2717: if (orientation >= 0) {
2718: const PetscInt dim = s->atlasDof[p];
2719: PetscInt i;
2721: if (mode == INSERT_VALUES) {
2722: for (i = 0; i < dim; ++i) array[i] = values ? values[i] : i;
2723: } else {
2724: for (i = 0; i < dim; ++i) array[i] += values[i];
2725: }
2726: } else {
2727: PetscInt offset = 0;
2728: PetscInt j = -1, field, i;
2730: for (field = 0; field < s->numFields; ++field) {
2731: const PetscInt dim = s->field[field]->atlasDof[p];
2733: for (i = dim - 1; i >= 0; --i) array[++j] = values ? values[i + offset] : i + offset;
2734: offset += dim;
2735: }
2736: }
2737: } else {
2738: if (orientation >= 0) {
2739: const PetscInt dim = s->atlasDof[p];
2740: PetscInt cInd = 0, i;
2741: const PetscInt *cDof;
2743: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2744: if (mode == INSERT_VALUES) {
2745: for (i = 0; i < dim; ++i) {
2746: if ((cInd < cDim) && (i == cDof[cInd])) {
2747: ++cInd;
2748: continue;
2749: }
2750: array[i] = values ? values[i] : i;
2751: }
2752: } else {
2753: for (i = 0; i < dim; ++i) {
2754: if ((cInd < cDim) && (i == cDof[cInd])) {
2755: ++cInd;
2756: continue;
2757: }
2758: array[i] += values[i];
2759: }
2760: }
2761: } else {
2762: const PetscInt *cDof;
2763: PetscInt offset = 0;
2764: PetscInt cOffset = 0;
2765: PetscInt j = 0, field;
2767: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2768: for (field = 0; field < s->numFields; ++field) {
2769: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2770: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2771: const PetscInt sDim = dim - tDim;
2772: PetscInt cInd = 0, i, k;
2774: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2775: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2776: ++cInd;
2777: continue;
2778: }
2779: array[j] = values ? values[k] : k;
2780: }
2781: offset += dim;
2782: cOffset += dim - tDim;
2783: }
2784: }
2785: }
2786: PetscFunctionReturn(PETSC_SUCCESS);
2787: }
2789: /*@
2790: PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs
2792: Not Collective
2794: Input Parameter:
2795: . s - The `PetscSection`
2797: Output Parameter:
2798: . hasConstraints - flag indicating that the section has constrained dofs
2800: Level: intermediate
2802: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2803: @*/
2804: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2805: {
2806: PetscFunctionBegin;
2808: PetscAssertPointer(hasConstraints, 2);
2809: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2810: PetscFunctionReturn(PETSC_SUCCESS);
2811: }
2813: /*@C
2814: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained for a given point
2816: Not Collective
2818: Input Parameters:
2819: + s - The `PetscSection`
2820: - point - The point
2822: Output Parameter:
2823: . indices - The constrained dofs
2825: Level: intermediate
2827: Fortran Notes:
2828: Use `PetscSectionGetConstraintIndicesF90()` and `PetscSectionRestoreConstraintIndicesF90()`
2830: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2831: @*/
2832: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt *indices[])
2833: {
2834: PetscFunctionBegin;
2836: if (s->bc) {
2837: PetscCall(VecIntGetValuesSection_Private(s->bcIndices, s->bc, point, indices));
2838: } else *indices = NULL;
2839: PetscFunctionReturn(PETSC_SUCCESS);
2840: }
2842: /*@
2843: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2845: Not Collective
2847: Input Parameters:
2848: + s - The `PetscSection`
2849: . point - The point
2850: - indices - The constrained dofs
2852: Level: intermediate
2854: Fortran Notes:
2855: Use `PetscSectionSetConstraintIndicesF90()`
2857: .seealso: [PetscSection](sec_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2858: @*/
2859: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2860: {
2861: PetscFunctionBegin;
2863: if (s->bc) {
2864: const PetscInt dof = s->atlasDof[point];
2865: const PetscInt cdof = s->bc->atlasDof[point];
2866: PetscInt d;
2868: if (indices)
2869: 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]);
2870: PetscCall(VecIntSetValuesSection_Private(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
2871: }
2872: PetscFunctionReturn(PETSC_SUCCESS);
2873: }
2875: /*@C
2876: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
2878: Not Collective
2880: Input Parameters:
2881: + s - The `PetscSection`
2882: . field - The field number
2883: - point - The point
2885: Output Parameter:
2886: . indices - The constrained dofs sorted in ascending order
2888: Level: intermediate
2890: Note:
2891: The indices array, which is provided by the caller, must have capacity to hold the number of constrained dofs, e.g., as returned by `PetscSectionGetConstraintDof()`.
2893: Fortran Notes:
2894: Use `PetscSectionGetFieldConstraintIndicesF90()` and `PetscSectionRestoreFieldConstraintIndicesF90()`
2896: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2897: @*/
2898: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt **indices)
2899: {
2900: PetscFunctionBegin;
2902: PetscAssertPointer(indices, 4);
2903: PetscSectionCheckValidField(field, s->numFields);
2904: PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
2905: PetscFunctionReturn(PETSC_SUCCESS);
2906: }
2908: /*@C
2909: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
2911: Not Collective
2913: Input Parameters:
2914: + s - The `PetscSection`
2915: . point - The point
2916: . field - The field number
2917: - indices - The constrained dofs
2919: Level: intermediate
2921: Fortran Notes:
2922: Use `PetscSectionSetFieldConstraintIndicesF90()`
2924: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2925: @*/
2926: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
2927: {
2928: PetscFunctionBegin;
2930: PetscSectionCheckValidField(field, s->numFields);
2931: PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
2932: PetscFunctionReturn(PETSC_SUCCESS);
2933: }
2935: /*@
2936: PetscSectionPermute - Reorder the section according to the input point permutation
2938: Collective
2940: Input Parameters:
2941: + section - The `PetscSection` object
2942: - permutation - The point permutation, old point p becomes new point perm[p]
2944: Output Parameter:
2945: . sectionNew - The permuted `PetscSection`
2947: Level: intermediate
2949: Note:
2950: The data and the access to the data via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` are both changed in `sectionNew`
2952: Compare to `PetscSectionSetPermutation()`
2954: .seealso: [PetscSection](sec_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
2955: @*/
2956: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
2957: {
2958: PetscSection s = section, sNew;
2959: const PetscInt *perm;
2960: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
2962: PetscFunctionBegin;
2965: PetscAssertPointer(sectionNew, 3);
2966: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
2967: PetscCall(PetscSectionGetNumFields(s, &numFields));
2968: if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
2969: for (f = 0; f < numFields; ++f) {
2970: const char *name;
2971: PetscInt numComp;
2973: PetscCall(PetscSectionGetFieldName(s, f, &name));
2974: PetscCall(PetscSectionSetFieldName(sNew, f, name));
2975: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2976: PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
2977: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2978: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2979: PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
2980: }
2981: }
2982: PetscCall(ISGetLocalSize(permutation, &numPoints));
2983: PetscCall(ISGetIndices(permutation, &perm));
2984: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2985: PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
2986: PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
2987: for (p = pStart; p < pEnd; ++p) {
2988: PetscInt dof, cdof;
2990: PetscCall(PetscSectionGetDof(s, p, &dof));
2991: PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
2992: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2993: if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
2994: for (f = 0; f < numFields; ++f) {
2995: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
2996: PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
2997: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
2998: if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
2999: }
3000: }
3001: PetscCall(PetscSectionSetUp(sNew));
3002: for (p = pStart; p < pEnd; ++p) {
3003: const PetscInt *cind;
3004: PetscInt cdof;
3006: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
3007: if (cdof) {
3008: PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
3009: PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
3010: }
3011: for (f = 0; f < numFields; ++f) {
3012: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
3013: if (cdof) {
3014: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
3015: PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
3016: }
3017: }
3018: }
3019: PetscCall(ISRestoreIndices(permutation, &perm));
3020: *sectionNew = sNew;
3021: PetscFunctionReturn(PETSC_SUCCESS);
3022: }
3024: /*@
3025: PetscSectionSetClosureIndex - Create an internal data structure to speed up closure queries.
3027: Collective
3029: Input Parameters:
3030: + section - The `PetscSection`
3031: . obj - A `PetscObject` which serves as the key for this index
3032: . clSection - `PetscSection` giving the size of the closure of each point
3033: - clPoints - `IS` giving the points in each closure
3035: Level: advanced
3037: Note:
3038: This function creates an internal map from each point to its closure. We compress out closure points with no dofs in this section.
3040: Developer Notes:
3041: The information provided here is completely opaque
3043: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
3044: @*/
3045: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
3046: {
3047: PetscFunctionBegin;
3051: if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
3052: section->clObj = obj;
3053: PetscCall(PetscObjectReference((PetscObject)clSection));
3054: PetscCall(PetscObjectReference((PetscObject)clPoints));
3055: PetscCall(PetscSectionDestroy(§ion->clSection));
3056: PetscCall(ISDestroy(§ion->clPoints));
3057: section->clSection = clSection;
3058: section->clPoints = clPoints;
3059: PetscFunctionReturn(PETSC_SUCCESS);
3060: }
3062: /*@
3063: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
3065: Collective
3067: Input Parameters:
3068: + section - The `PetscSection`
3069: - obj - A `PetscObject` which serves as the key for this index
3071: Output Parameters:
3072: + clSection - `PetscSection` giving the size of the closure of each point
3073: - clPoints - `IS` giving the points in each closure
3075: Level: advanced
3077: .seealso: [PetscSection](sec_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3078: @*/
3079: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
3080: {
3081: PetscFunctionBegin;
3082: if (section->clObj == obj) {
3083: if (clSection) *clSection = section->clSection;
3084: if (clPoints) *clPoints = section->clPoints;
3085: } else {
3086: if (clSection) *clSection = NULL;
3087: if (clPoints) *clPoints = NULL;
3088: }
3089: PetscFunctionReturn(PETSC_SUCCESS);
3090: }
3092: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
3093: {
3094: PetscInt i;
3095: khiter_t iter;
3096: int new_entry;
3097: PetscSectionClosurePermKey key = {depth, clSize};
3098: PetscSectionClosurePermVal *val;
3100: PetscFunctionBegin;
3101: if (section->clObj != obj) {
3102: PetscCall(PetscSectionDestroy(§ion->clSection));
3103: PetscCall(ISDestroy(§ion->clPoints));
3104: }
3105: section->clObj = obj;
3106: if (!section->clHash) PetscCall(PetscClPermCreate(§ion->clHash));
3107: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
3108: val = &kh_val(section->clHash, iter);
3109: if (!new_entry) {
3110: PetscCall(PetscFree(val->perm));
3111: PetscCall(PetscFree(val->invPerm));
3112: }
3113: if (mode == PETSC_COPY_VALUES) {
3114: PetscCall(PetscMalloc1(clSize, &val->perm));
3115: PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
3116: } else if (mode == PETSC_OWN_POINTER) {
3117: val->perm = clPerm;
3118: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
3119: PetscCall(PetscMalloc1(clSize, &val->invPerm));
3120: for (i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
3121: PetscFunctionReturn(PETSC_SUCCESS);
3122: }
3124: /*@
3125: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3127: Not Collective
3129: Input Parameters:
3130: + section - The `PetscSection`
3131: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3132: . depth - Depth of points on which to apply the given permutation
3133: - perm - Permutation of the cell dof closure
3135: Level: intermediate
3137: Notes:
3138: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
3139: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
3140: each topology and degree.
3142: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
3143: exotic/enriched spaces on mixed topology meshes.
3145: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
3146: @*/
3147: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
3148: {
3149: const PetscInt *clPerm = NULL;
3150: PetscInt clSize = 0;
3152: PetscFunctionBegin;
3153: if (perm) {
3154: PetscCall(ISGetLocalSize(perm, &clSize));
3155: PetscCall(ISGetIndices(perm, &clPerm));
3156: }
3157: PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
3158: if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
3159: PetscFunctionReturn(PETSC_SUCCESS);
3160: }
3162: static PetscErrorCode PetscSectionGetClosurePermutation_Private(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3163: {
3164: PetscFunctionBegin;
3165: if (section->clObj == obj) {
3166: PetscSectionClosurePermKey k = {depth, size};
3167: PetscSectionClosurePermVal v;
3169: PetscCall(PetscClPermGet(section->clHash, k, &v));
3170: if (perm) *perm = v.perm;
3171: } else {
3172: if (perm) *perm = NULL;
3173: }
3174: PetscFunctionReturn(PETSC_SUCCESS);
3175: }
3177: /*@
3178: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3180: Not Collective
3182: Input Parameters:
3183: + section - The `PetscSection`
3184: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
3185: . depth - Depth stratum on which to obtain closure permutation
3186: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3188: Output Parameter:
3189: . perm - The dof closure permutation
3191: Level: intermediate
3193: Note:
3194: The user must destroy the `IS` that is returned.
3196: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3197: @*/
3198: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3199: {
3200: const PetscInt *clPerm = NULL;
3202: PetscFunctionBegin;
3203: PetscCall(PetscSectionGetClosurePermutation_Private(section, obj, depth, clSize, &clPerm));
3204: PetscCheck(clPerm, PetscObjectComm((PetscObject)obj), PETSC_ERR_ARG_WRONG, "There is no closure permutation associated with this object for depth %" PetscInt_FMT " of size %" PetscInt_FMT, depth, clSize);
3205: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3206: PetscFunctionReturn(PETSC_SUCCESS);
3207: }
3209: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3210: {
3211: PetscFunctionBegin;
3212: if (section->clObj == obj && section->clHash) {
3213: PetscSectionClosurePermKey k = {depth, size};
3214: PetscSectionClosurePermVal v;
3215: PetscCall(PetscClPermGet(section->clHash, k, &v));
3216: if (perm) *perm = v.invPerm;
3217: } else {
3218: if (perm) *perm = NULL;
3219: }
3220: PetscFunctionReturn(PETSC_SUCCESS);
3221: }
3223: /*@
3224: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
3226: Not Collective
3228: Input Parameters:
3229: + section - The `PetscSection`
3230: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3231: . depth - Depth stratum on which to obtain closure permutation
3232: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3234: Output Parameter:
3235: . perm - The dof closure permutation
3237: Level: intermediate
3239: Note:
3240: The user must destroy the `IS` that is returned.
3242: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3243: @*/
3244: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3245: {
3246: const PetscInt *clPerm = NULL;
3248: PetscFunctionBegin;
3249: PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3250: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3251: PetscFunctionReturn(PETSC_SUCCESS);
3252: }
3254: /*@
3255: PetscSectionGetField - Get the `PetscSection` associated with a single field
3257: Input Parameters:
3258: + s - The `PetscSection`
3259: - field - The field number
3261: Output Parameter:
3262: . subs - The `PetscSection` for the given field, note the chart of `subs` is not set
3264: Level: intermediate
3266: Note:
3267: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
3269: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3270: @*/
3271: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3272: {
3273: PetscFunctionBegin;
3275: PetscAssertPointer(subs, 3);
3276: PetscSectionCheckValidField(field, s->numFields);
3277: *subs = s->field[field];
3278: PetscFunctionReturn(PETSC_SUCCESS);
3279: }
3281: PetscClassId PETSC_SECTION_SYM_CLASSID;
3282: PetscFunctionList PetscSectionSymList = NULL;
3284: /*@
3285: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
3287: Collective
3289: Input Parameter:
3290: . comm - the MPI communicator
3292: Output Parameter:
3293: . sym - pointer to the new set of symmetries
3295: Level: developer
3297: .seealso: [PetscSection](sec_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3298: @*/
3299: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3300: {
3301: PetscFunctionBegin;
3302: PetscAssertPointer(sym, 2);
3303: PetscCall(ISInitializePackage());
3305: PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3306: PetscFunctionReturn(PETSC_SUCCESS);
3307: }
3309: /*@
3310: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
3312: Collective
3314: Input Parameters:
3315: + sym - The section symmetry object
3316: - method - The name of the section symmetry type
3318: Level: developer
3320: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3321: @*/
3322: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3323: {
3324: PetscErrorCode (*r)(PetscSectionSym);
3325: PetscBool match;
3327: PetscFunctionBegin;
3329: PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3330: if (match) PetscFunctionReturn(PETSC_SUCCESS);
3332: PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3333: PetscCheck(r, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3334: PetscTryTypeMethod(sym, destroy);
3335: sym->ops->destroy = NULL;
3337: PetscCall((*r)(sym));
3338: PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3339: PetscFunctionReturn(PETSC_SUCCESS);
3340: }
3342: /*@
3343: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
3345: Not Collective
3347: Input Parameter:
3348: . sym - The section symmetry
3350: Output Parameter:
3351: . type - The index set type name
3353: Level: developer
3355: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3356: @*/
3357: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3358: {
3359: PetscFunctionBegin;
3361: PetscAssertPointer(type, 2);
3362: *type = ((PetscObject)sym)->type_name;
3363: PetscFunctionReturn(PETSC_SUCCESS);
3364: }
3366: /*@C
3367: PetscSectionSymRegister - Registers a new section symmetry implementation
3369: Not Collective, No Fortran Support
3371: Input Parameters:
3372: + sname - The name of a new user-defined creation routine
3373: - function - The creation routine itself
3375: Level: developer
3377: Notes:
3378: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
3380: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3381: @*/
3382: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3383: {
3384: PetscFunctionBegin;
3385: PetscCall(ISInitializePackage());
3386: PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3387: PetscFunctionReturn(PETSC_SUCCESS);
3388: }
3390: /*@
3391: PetscSectionSymDestroy - Destroys a section symmetry.
3393: Collective
3395: Input Parameter:
3396: . sym - the section symmetry
3398: Level: developer
3400: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`
3401: @*/
3402: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3403: {
3404: SymWorkLink link, next;
3406: PetscFunctionBegin;
3407: if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3409: if (--((PetscObject)*sym)->refct > 0) {
3410: *sym = NULL;
3411: PetscFunctionReturn(PETSC_SUCCESS);
3412: }
3413: PetscTryTypeMethod(*sym, destroy);
3414: PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3415: for (link = (*sym)->workin; link; link = next) {
3416: PetscInt **perms = (PetscInt **)link->perms;
3417: PetscScalar **rots = (PetscScalar **)link->rots;
3418: PetscCall(PetscFree2(perms, rots));
3419: next = link->next;
3420: PetscCall(PetscFree(link));
3421: }
3422: (*sym)->workin = NULL;
3423: PetscCall(PetscHeaderDestroy(sym));
3424: PetscFunctionReturn(PETSC_SUCCESS);
3425: }
3427: /*@
3428: PetscSectionSymView - Displays a section symmetry
3430: Collective
3432: Input Parameters:
3433: + sym - the index set
3434: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3436: Level: developer
3438: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3439: @*/
3440: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3441: {
3442: PetscFunctionBegin;
3444: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3446: PetscCheckSameComm(sym, 1, viewer, 2);
3447: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3448: PetscTryTypeMethod(sym, view, viewer);
3449: PetscFunctionReturn(PETSC_SUCCESS);
3450: }
3452: /*@
3453: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3455: Collective
3457: Input Parameters:
3458: + section - the section describing data layout
3459: - sym - the symmetry describing the affect of orientation on the access of the data
3461: Level: developer
3463: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3464: @*/
3465: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3466: {
3467: PetscFunctionBegin;
3469: PetscCall(PetscSectionSymDestroy(§ion->sym));
3470: if (sym) {
3472: PetscCheckSameComm(section, 1, sym, 2);
3473: PetscCall(PetscObjectReference((PetscObject)sym));
3474: }
3475: section->sym = sym;
3476: PetscFunctionReturn(PETSC_SUCCESS);
3477: }
3479: /*@
3480: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3482: Not Collective
3484: Input Parameter:
3485: . section - the section describing data layout
3487: Output Parameter:
3488: . sym - the symmetry describing the affect of orientation on the access of the data, provided previously by `PetscSectionSetSym()`
3490: Level: developer
3492: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3493: @*/
3494: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3495: {
3496: PetscFunctionBegin;
3498: *sym = section->sym;
3499: PetscFunctionReturn(PETSC_SUCCESS);
3500: }
3502: /*@
3503: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3505: Collective
3507: Input Parameters:
3508: + section - the section describing data layout
3509: . field - the field number
3510: - sym - the symmetry describing the affect of orientation on the access of the data
3512: Level: developer
3514: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3515: @*/
3516: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3517: {
3518: PetscFunctionBegin;
3520: PetscSectionCheckValidField(field, section->numFields);
3521: PetscCall(PetscSectionSetSym(section->field[field], sym));
3522: PetscFunctionReturn(PETSC_SUCCESS);
3523: }
3525: /*@
3526: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3528: Collective
3530: Input Parameters:
3531: + section - the section describing data layout
3532: - field - the field number
3534: Output Parameter:
3535: . sym - the symmetry describing the affect of orientation on the access of the data
3537: Level: developer
3539: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3540: @*/
3541: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3542: {
3543: PetscFunctionBegin;
3545: PetscSectionCheckValidField(field, section->numFields);
3546: *sym = section->field[field]->sym;
3547: PetscFunctionReturn(PETSC_SUCCESS);
3548: }
3550: /*@C
3551: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3553: Not Collective
3555: Input Parameters:
3556: + section - the section
3557: . numPoints - the number of points
3558: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3559: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3560: context, see `DMPlexGetConeOrientation()`).
3562: Output Parameters:
3563: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3564: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3565: identity).
3567: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3568: .vb
3569: const PetscInt **perms;
3570: const PetscScalar **rots;
3571: PetscInt lOffset;
3573: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3574: for (i = 0, lOffset = 0; i < numPoints; i++) {
3575: PetscInt point = points[2*i], dof, sOffset;
3576: const PetscInt *perm = perms ? perms[i] : NULL;
3577: const PetscScalar *rot = rots ? rots[i] : NULL;
3579: PetscSectionGetDof(section,point,&dof);
3580: PetscSectionGetOffset(section,point,&sOffset);
3582: if (perm) {for (j = 0; j < dof; j++) {lArray[lOffset + perm[j]] = sArray[sOffset + j];}}
3583: else {for (j = 0; j < dof; j++) {lArray[lOffset + j ] = sArray[sOffset + j];}}
3584: if (rot) {for (j = 0; j < dof; j++) {lArray[lOffset + j ] *= rot[j]; }}
3585: lOffset += dof;
3586: }
3587: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3588: .ve
3590: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3591: .vb
3592: const PetscInt **perms;
3593: const PetscScalar **rots;
3594: PetscInt lOffset;
3596: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3597: for (i = 0, lOffset = 0; i < numPoints; i++) {
3598: PetscInt point = points[2*i], dof, sOffset;
3599: const PetscInt *perm = perms ? perms[i] : NULL;
3600: const PetscScalar *rot = rots ? rots[i] : NULL;
3602: PetscSectionGetDof(section,point,&dof);
3603: PetscSectionGetOffset(section,point,&sOff);
3605: if (perm) {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.);}}
3606: else {for (j = 0; j < dof; j++) {sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.);}}
3607: offset += dof;
3608: }
3609: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3610: .ve
3612: Level: developer
3614: Notes:
3615: `PetscSectionSetSym()` must have been previously called to provide the symmetries to the `PetscSection`
3617: Use `PetscSectionRestorePointSyms()` when finished with the data
3619: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3620: @*/
3621: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3622: {
3623: PetscSectionSym sym;
3625: PetscFunctionBegin;
3627: if (numPoints) PetscAssertPointer(points, 3);
3628: if (perms) *perms = NULL;
3629: if (rots) *rots = NULL;
3630: sym = section->sym;
3631: if (sym && (perms || rots)) {
3632: SymWorkLink link;
3634: if (sym->workin) {
3635: link = sym->workin;
3636: sym->workin = sym->workin->next;
3637: } else {
3638: PetscCall(PetscNew(&link));
3639: }
3640: if (numPoints > link->numPoints) {
3641: PetscInt **perms = (PetscInt **)link->perms;
3642: PetscScalar **rots = (PetscScalar **)link->rots;
3643: PetscCall(PetscFree2(perms, rots));
3644: PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3645: link->numPoints = numPoints;
3646: }
3647: link->next = sym->workout;
3648: sym->workout = link;
3649: PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3650: PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3651: PetscUseTypeMethod(sym, getpoints, section, numPoints, points, link->perms, link->rots);
3652: if (perms) *perms = link->perms;
3653: if (rots) *rots = link->rots;
3654: }
3655: PetscFunctionReturn(PETSC_SUCCESS);
3656: }
3658: /*@C
3659: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3661: Not Collective
3663: Input Parameters:
3664: + section - the section
3665: . numPoints - the number of points
3666: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3667: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3668: context, see `DMPlexGetConeOrientation()`).
3669: . perms - The permutations for the given orientations: set to `NULL` at conclusion
3670: - rots - The field rotations symmetries for the given orientations: set to `NULL` at conclusion
3672: Level: developer
3674: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3675: @*/
3676: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3677: {
3678: PetscSectionSym sym;
3680: PetscFunctionBegin;
3682: sym = section->sym;
3683: if (sym && (perms || rots)) {
3684: SymWorkLink *p, link;
3686: for (p = &sym->workout; (link = *p); p = &link->next) {
3687: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3688: *p = link->next;
3689: link->next = sym->workin;
3690: sym->workin = link;
3691: if (perms) *perms = NULL;
3692: if (rots) *rots = NULL;
3693: PetscFunctionReturn(PETSC_SUCCESS);
3694: }
3695: }
3696: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3697: }
3698: PetscFunctionReturn(PETSC_SUCCESS);
3699: }
3701: /*@C
3702: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3704: Not Collective
3706: Input Parameters:
3707: + section - the section
3708: . field - the field of the section
3709: . numPoints - the number of points
3710: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3711: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3712: context, see `DMPlexGetConeOrientation()`).
3714: Output Parameters:
3715: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3716: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3717: identity).
3719: Level: developer
3721: Notes:
3722: `PetscSectionSetFieldSym()` must have been previously called to provide the symmetries to the `PetscSection`
3724: Use `PetscSectionRestoreFieldPointSyms()` when finished with the data
3726: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3727: @*/
3728: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3729: {
3730: PetscFunctionBegin;
3732: 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);
3733: PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3734: PetscFunctionReturn(PETSC_SUCCESS);
3735: }
3737: /*@C
3738: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3740: Not Collective
3742: Input Parameters:
3743: + section - the section
3744: . field - the field number
3745: . numPoints - the number of points
3746: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3747: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3748: context, see `DMPlexGetConeOrientation()`).
3749: . perms - The permutations for the given orientations: set to NULL at conclusion
3750: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3752: Level: developer
3754: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3755: @*/
3756: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3757: {
3758: PetscFunctionBegin;
3760: 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);
3761: PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3762: PetscFunctionReturn(PETSC_SUCCESS);
3763: }
3765: /*@
3766: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3768: Not Collective
3770: Input Parameter:
3771: . sym - the `PetscSectionSym`
3773: Output Parameter:
3774: . nsym - the equivalent symmetries
3776: Level: developer
3778: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3779: @*/
3780: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3781: {
3782: PetscFunctionBegin;
3785: PetscTryTypeMethod(sym, copy, nsym);
3786: PetscFunctionReturn(PETSC_SUCCESS);
3787: }
3789: /*@
3790: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3792: Collective
3794: Input Parameters:
3795: + sym - the `PetscSectionSym`
3796: - migrationSF - the distribution map from roots to leaves
3798: Output Parameter:
3799: . dsym - the redistributed symmetries
3801: Level: developer
3803: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3804: @*/
3805: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3806: {
3807: PetscFunctionBegin;
3810: PetscAssertPointer(dsym, 3);
3811: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3812: PetscFunctionReturn(PETSC_SUCCESS);
3813: }
3815: /*@
3816: PetscSectionGetUseFieldOffsets - Get the flag indicating if field offsets are used directly in a global section, rather than just the point offset
3818: Not Collective
3820: Input Parameter:
3821: . s - the global `PetscSection`
3823: Output Parameter:
3824: . flg - the flag
3826: Level: developer
3828: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3829: @*/
3830: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3831: {
3832: PetscFunctionBegin;
3834: *flg = s->useFieldOff;
3835: PetscFunctionReturn(PETSC_SUCCESS);
3836: }
3838: /*@
3839: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3841: Not Collective
3843: Input Parameters:
3844: + s - the global `PetscSection`
3845: - flg - the flag
3847: Level: developer
3849: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3850: @*/
3851: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3852: {
3853: PetscFunctionBegin;
3855: s->useFieldOff = flg;
3856: PetscFunctionReturn(PETSC_SUCCESS);
3857: }
3859: #define PetscSectionExpandPoints_Loop(TYPE) \
3860: do { \
3861: PetscInt i, n, o0, o1, size; \
3862: TYPE *a0 = (TYPE *)origArray, *a1; \
3863: PetscCall(PetscSectionGetStorageSize(s, &size)); \
3864: PetscCall(PetscMalloc1(size, &a1)); \
3865: for (i = 0; i < npoints; i++) { \
3866: PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
3867: PetscCall(PetscSectionGetOffset(s, i, &o1)); \
3868: PetscCall(PetscSectionGetDof(s, i, &n)); \
3869: PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n *unitsize)); \
3870: } \
3871: *newArray = (void *)a1; \
3872: } while (0)
3874: /*@
3875: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3877: Not Collective
3879: Input Parameters:
3880: + origSection - the `PetscSection` describing the layout of the array
3881: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
3882: . origArray - the array; its size must be equal to the storage size of `origSection`
3883: - points - `IS` with points to extract; its indices must lie in the chart of `origSection`
3885: Output Parameters:
3886: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
3887: - newArray - the array of the extracted DOFs; its size is the storage size of `newSection`
3889: Level: developer
3891: .seealso: [PetscSection](sec_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
3892: @*/
3893: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
3894: {
3895: PetscSection s;
3896: const PetscInt *points_;
3897: PetscInt i, n, npoints, pStart, pEnd;
3898: PetscMPIInt unitsize;
3900: PetscFunctionBegin;
3902: PetscAssertPointer(origArray, 3);
3904: if (newSection) PetscAssertPointer(newSection, 5);
3905: if (newArray) PetscAssertPointer(newArray, 6);
3906: PetscCallMPI(MPI_Type_size(dataType, &unitsize));
3907: PetscCall(ISGetLocalSize(points, &npoints));
3908: PetscCall(ISGetIndices(points, &points_));
3909: PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
3910: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
3911: PetscCall(PetscSectionSetChart(s, 0, npoints));
3912: for (i = 0; i < npoints; i++) {
3913: 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);
3914: PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
3915: PetscCall(PetscSectionSetDof(s, i, n));
3916: }
3917: PetscCall(PetscSectionSetUp(s));
3918: if (newArray) {
3919: if (dataType == MPIU_INT) {
3920: PetscSectionExpandPoints_Loop(PetscInt);
3921: } else if (dataType == MPIU_SCALAR) {
3922: PetscSectionExpandPoints_Loop(PetscScalar);
3923: } else if (dataType == MPIU_REAL) {
3924: PetscSectionExpandPoints_Loop(PetscReal);
3925: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
3926: }
3927: if (newSection) {
3928: *newSection = s;
3929: } else {
3930: PetscCall(PetscSectionDestroy(&s));
3931: }
3932: PetscCall(ISRestoreIndices(points, &points_));
3933: PetscFunctionReturn(PETSC_SUCCESS);
3934: }