Actual source code: section.c
1: /*
2: This file contains routines for basic section object implementation.
3: */
5: #include <petsc/private/sectionimpl.h>
6: #include <petscsf.h>
8: PetscClassId PETSC_SECTION_CLASSID;
10: /*@
11: PetscSectionCreate - Allocates a `PetscSection` and sets the map contents to the default.
13: Collective
15: Input Parameters:
16: + comm - the MPI communicator
17: - s - pointer to the section
19: Level: beginner
21: Notes:
22: Typical calling sequence
23: .vb
24: PetscSectionCreate(MPI_Comm,PetscSection *);!
25: PetscSectionSetNumFields(PetscSection, numFields);
26: PetscSectionSetChart(PetscSection,low,high);
27: PetscSectionSetDof(PetscSection,point,numdof);
28: PetscSectionSetUp(PetscSection);
29: PetscSectionGetOffset(PetscSection,point,PetscInt *);
30: PetscSectionDestroy(PetscSection);
31: .ve
33: The `PetscSection` object and methods are intended to be used in the PETSc `Vec` and `Mat` implementations. The indices returned by the `PetscSection` are appropriate for the kind of `Vec` it is associated with. For example, if the vector being indexed is a local vector, we call the section a local section. If the section indexes a global vector, we call it a global section. For parallel vectors, like global vectors, we use negative indices to indicate dofs owned by other processes.
35: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetChart()`, `PetscSectionDestroy()`, `PetscSectionCreateGlobalSection()`
36: @*/
37: PetscErrorCode PetscSectionCreate(MPI_Comm comm, PetscSection *s)
38: {
39: PetscFunctionBegin;
40: PetscAssertPointer(s, 2);
41: PetscCall(ISInitializePackage());
43: PetscCall(PetscHeaderCreate(*s, PETSC_SECTION_CLASSID, "PetscSection", "Section", "IS", comm, PetscSectionDestroy, PetscSectionView));
44: (*s)->pStart = -1;
45: (*s)->pEnd = -1;
46: (*s)->perm = NULL;
47: (*s)->pointMajor = PETSC_TRUE;
48: (*s)->includesConstraints = PETSC_TRUE;
49: (*s)->atlasDof = NULL;
50: (*s)->atlasOff = NULL;
51: (*s)->bc = NULL;
52: (*s)->bcIndices = NULL;
53: (*s)->setup = PETSC_FALSE;
54: (*s)->numFields = 0;
55: (*s)->fieldNames = NULL;
56: (*s)->field = NULL;
57: (*s)->useFieldOff = PETSC_FALSE;
58: (*s)->compNames = NULL;
59: (*s)->clObj = NULL;
60: (*s)->clHash = NULL;
61: (*s)->clSection = NULL;
62: (*s)->clPoints = NULL;
63: PetscCall(PetscSectionInvalidateMaxDof_Internal(*s));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*@
68: PetscSectionCopy - Creates a shallow (if possible) copy of the `PetscSection`
70: Collective
72: Input Parameter:
73: . section - the `PetscSection`
75: Output Parameter:
76: . newSection - the copy
78: Level: intermediate
80: Developer Notes:
81: What exactly does shallow mean in this context?
83: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
84: @*/
85: PetscErrorCode PetscSectionCopy(PetscSection section, PetscSection newSection)
86: {
87: PetscFunctionBegin;
90: PetscCall(PetscSectionCopy_Internal(section, newSection, NULL));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: PetscErrorCode PetscSectionCopy_Internal(PetscSection section, PetscSection newSection, PetscBT constrained_dofs)
95: {
96: PetscSectionSym sym;
97: IS perm;
98: PetscInt numFields, f, c, pStart, pEnd, p;
100: PetscFunctionBegin;
103: PetscCall(PetscSectionReset(newSection));
104: PetscCall(PetscSectionGetNumFields(section, &numFields));
105: if (numFields) PetscCall(PetscSectionSetNumFields(newSection, numFields));
106: for (f = 0; f < numFields; ++f) {
107: const char *fieldName = NULL, *compName = NULL;
108: PetscInt numComp = 0;
110: PetscCall(PetscSectionGetFieldName(section, f, &fieldName));
111: PetscCall(PetscSectionSetFieldName(newSection, f, fieldName));
112: PetscCall(PetscSectionGetFieldComponents(section, f, &numComp));
113: PetscCall(PetscSectionSetFieldComponents(newSection, f, numComp));
114: for (c = 0; c < numComp; ++c) {
115: PetscCall(PetscSectionGetComponentName(section, f, c, &compName));
116: PetscCall(PetscSectionSetComponentName(newSection, f, c, compName));
117: }
118: PetscCall(PetscSectionGetFieldSym(section, f, &sym));
119: PetscCall(PetscSectionSetFieldSym(newSection, f, sym));
120: }
121: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
122: PetscCall(PetscSectionSetChart(newSection, pStart, pEnd));
123: PetscCall(PetscSectionGetPermutation(section, &perm));
124: PetscCall(PetscSectionSetPermutation(newSection, perm));
125: PetscCall(PetscSectionGetSym(section, &sym));
126: PetscCall(PetscSectionSetSym(newSection, sym));
127: for (p = pStart; p < pEnd; ++p) {
128: PetscInt dof, cdof, fcdof = 0;
129: PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));
131: PetscCall(PetscSectionGetDof(section, p, &dof));
132: PetscCall(PetscSectionSetDof(newSection, p, dof));
133: if (force_constrained) cdof = dof;
134: else PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
135: if (cdof) PetscCall(PetscSectionSetConstraintDof(newSection, p, cdof));
136: for (f = 0; f < numFields; ++f) {
137: PetscCall(PetscSectionGetFieldDof(section, p, f, &dof));
138: PetscCall(PetscSectionSetFieldDof(newSection, p, f, dof));
139: if (cdof) {
140: if (force_constrained) fcdof = dof;
141: else PetscCall(PetscSectionGetFieldConstraintDof(section, p, f, &fcdof));
142: if (fcdof) PetscCall(PetscSectionSetFieldConstraintDof(newSection, p, f, fcdof));
143: }
144: }
145: }
146: PetscCall(PetscSectionSetUp(newSection));
147: for (p = pStart; p < pEnd; ++p) {
148: PetscInt off, cdof, fcdof = 0;
149: const PetscInt *cInd;
150: PetscBool force_constrained = (PetscBool)(constrained_dofs && PetscBTLookup(constrained_dofs, p - pStart));
152: /* Must set offsets in case they do not agree with the prefix sums */
153: PetscCall(PetscSectionGetOffset(section, p, &off));
154: PetscCall(PetscSectionSetOffset(newSection, p, off));
155: PetscCall(PetscSectionGetConstraintDof(newSection, p, &cdof));
156: if (cdof) {
157: if (force_constrained) cInd = NULL;
158: else PetscCall(PetscSectionGetConstraintIndices(section, p, &cInd));
159: PetscCall(PetscSectionSetConstraintIndices(newSection, p, cInd));
160: for (f = 0; f < numFields; ++f) {
161: PetscCall(PetscSectionGetFieldOffset(section, p, f, &off));
162: PetscCall(PetscSectionSetFieldOffset(newSection, p, f, off));
163: PetscCall(PetscSectionGetFieldConstraintDof(newSection, p, f, &fcdof));
164: if (fcdof) {
165: if (force_constrained) cInd = NULL;
166: else PetscCall(PetscSectionGetFieldConstraintIndices(section, p, f, &cInd));
167: PetscCall(PetscSectionSetFieldConstraintIndices(newSection, p, f, cInd));
168: }
169: }
170: }
171: }
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: /*@
176: PetscSectionClone - Creates a shallow (if possible) copy of the `PetscSection`
178: Collective
180: Input Parameter:
181: . section - the `PetscSection`
183: Output Parameter:
184: . newSection - the copy
186: Level: beginner
188: Developer Notes:
189: With standard PETSc terminology this should be called `PetscSectionDuplicate()`
191: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionCopy()`
192: @*/
193: PetscErrorCode PetscSectionClone(PetscSection section, PetscSection *newSection)
194: {
195: PetscFunctionBegin;
197: PetscAssertPointer(newSection, 2);
198: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), newSection));
199: PetscCall(PetscSectionCopy(section, *newSection));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: PetscSectionSetFromOptions - sets parameters in a `PetscSection` from the options database
206: Collective
208: Input Parameter:
209: . s - the `PetscSection`
211: Options Database Key:
212: . -petscsection_point_major - `PETSC_TRUE` for point-major order
214: Level: intermediate
216: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`
217: @*/
218: PetscErrorCode PetscSectionSetFromOptions(PetscSection s)
219: {
220: PetscFunctionBegin;
222: PetscObjectOptionsBegin((PetscObject)s);
223: PetscCall(PetscOptionsBool("-petscsection_point_major", "The for ordering, either point major or field major", "PetscSectionSetPointMajor", s->pointMajor, &s->pointMajor, NULL));
224: /* process any options handlers added with PetscObjectAddOptionsHandler() */
225: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)s, PetscOptionsObject));
226: PetscOptionsEnd();
227: PetscCall(PetscObjectViewFromOptions((PetscObject)s, NULL, "-petscsection_view"));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: PetscSectionCompare - Compares two sections
234: Collective
236: Input Parameters:
237: + s1 - the first `PetscSection`
238: - s2 - the second `PetscSection`
240: Output Parameter:
241: . congruent - `PETSC_TRUE` if the two sections are congruent, `PETSC_FALSE` otherwise
243: Level: intermediate
245: Note:
246: Field names are disregarded.
248: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionCopy()`, `PetscSectionClone()`
249: @*/
250: PetscErrorCode PetscSectionCompare(PetscSection s1, PetscSection s2, PetscBool *congruent)
251: {
252: PetscInt pStart, pEnd, nfields, ncdof, nfcdof, p, f, n1, n2;
253: const PetscInt *idx1, *idx2;
254: IS perm1, perm2;
255: PetscBool flg;
256: PetscMPIInt mflg;
258: PetscFunctionBegin;
261: PetscAssertPointer(congruent, 3);
262: flg = PETSC_FALSE;
264: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)s1), PetscObjectComm((PetscObject)s2), &mflg));
265: if (mflg != MPI_CONGRUENT && mflg != MPI_IDENT) {
266: *congruent = PETSC_FALSE;
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: PetscCall(PetscSectionGetChart(s1, &pStart, &pEnd));
271: PetscCall(PetscSectionGetChart(s2, &n1, &n2));
272: if (pStart != n1 || pEnd != n2) goto not_congruent;
274: PetscCall(PetscSectionGetPermutation(s1, &perm1));
275: PetscCall(PetscSectionGetPermutation(s2, &perm2));
276: if (perm1 && perm2) {
277: PetscCall(ISEqual(perm1, perm2, congruent));
278: if (!*congruent) goto not_congruent;
279: } else if (perm1 != perm2) goto not_congruent;
281: for (p = pStart; p < pEnd; ++p) {
282: PetscCall(PetscSectionGetOffset(s1, p, &n1));
283: PetscCall(PetscSectionGetOffset(s2, p, &n2));
284: if (n1 != n2) goto not_congruent;
286: PetscCall(PetscSectionGetDof(s1, p, &n1));
287: PetscCall(PetscSectionGetDof(s2, p, &n2));
288: if (n1 != n2) goto not_congruent;
290: PetscCall(PetscSectionGetConstraintDof(s1, p, &ncdof));
291: PetscCall(PetscSectionGetConstraintDof(s2, p, &n2));
292: if (ncdof != n2) goto not_congruent;
294: PetscCall(PetscSectionGetConstraintIndices(s1, p, &idx1));
295: PetscCall(PetscSectionGetConstraintIndices(s2, p, &idx2));
296: PetscCall(PetscArraycmp(idx1, idx2, ncdof, congruent));
297: if (!*congruent) goto not_congruent;
298: }
300: PetscCall(PetscSectionGetNumFields(s1, &nfields));
301: PetscCall(PetscSectionGetNumFields(s2, &n2));
302: if (nfields != n2) goto not_congruent;
304: for (f = 0; f < nfields; ++f) {
305: PetscCall(PetscSectionGetFieldComponents(s1, f, &n1));
306: PetscCall(PetscSectionGetFieldComponents(s2, f, &n2));
307: if (n1 != n2) goto not_congruent;
309: for (p = pStart; p < pEnd; ++p) {
310: PetscCall(PetscSectionGetFieldOffset(s1, p, f, &n1));
311: PetscCall(PetscSectionGetFieldOffset(s2, p, f, &n2));
312: if (n1 != n2) goto not_congruent;
314: PetscCall(PetscSectionGetFieldDof(s1, p, f, &n1));
315: PetscCall(PetscSectionGetFieldDof(s2, p, f, &n2));
316: if (n1 != n2) goto not_congruent;
318: PetscCall(PetscSectionGetFieldConstraintDof(s1, p, f, &nfcdof));
319: PetscCall(PetscSectionGetFieldConstraintDof(s2, p, f, &n2));
320: if (nfcdof != n2) goto not_congruent;
322: PetscCall(PetscSectionGetFieldConstraintIndices(s1, p, f, &idx1));
323: PetscCall(PetscSectionGetFieldConstraintIndices(s2, p, f, &idx2));
324: PetscCall(PetscArraycmp(idx1, idx2, nfcdof, congruent));
325: if (!*congruent) goto not_congruent;
326: }
327: }
329: flg = PETSC_TRUE;
330: not_congruent:
331: PetscCallMPI(MPIU_Allreduce(&flg, congruent, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)s1)));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: /*@
336: PetscSectionGetNumFields - Returns the number of fields in a `PetscSection`, or 0 if no fields were defined.
338: Not Collective
340: Input Parameter:
341: . s - the `PetscSection`
343: Output Parameter:
344: . numFields - the number of fields defined, or 0 if none were defined
346: Level: intermediate
348: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetNumFields()`
349: @*/
350: PetscErrorCode PetscSectionGetNumFields(PetscSection s, PetscInt *numFields)
351: {
352: PetscFunctionBegin;
354: PetscAssertPointer(numFields, 2);
355: *numFields = s->numFields;
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: /*@
360: PetscSectionSetNumFields - Sets the number of fields in a `PetscSection`
362: Not Collective
364: Input Parameters:
365: + s - the `PetscSection`
366: - numFields - the number of fields
368: Level: intermediate
370: Notes:
371: Calling this destroys all the information in the `PetscSection` including the chart.
373: You must call `PetscSectionSetChart()` after calling this.
375: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetNumFields()`, `PetscSectionSetChart()`, `PetscSectionReset()`
376: @*/
377: PetscErrorCode PetscSectionSetNumFields(PetscSection s, PetscInt numFields)
378: {
379: PetscFunctionBegin;
381: PetscCheck(numFields > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The number of fields %" PetscInt_FMT " must be positive", numFields);
382: PetscCall(PetscSectionReset(s));
384: s->numFields = numFields;
385: PetscCall(PetscMalloc1(s->numFields, &s->numFieldComponents));
386: PetscCall(PetscMalloc1(s->numFields, &s->fieldNames));
387: PetscCall(PetscMalloc1(s->numFields, &s->compNames));
388: PetscCall(PetscMalloc1(s->numFields, &s->field));
389: for (PetscInt f = 0; f < s->numFields; ++f) {
390: char name[64];
392: s->numFieldComponents[f] = 1;
394: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &s->field[f]));
395: PetscCall(PetscSNPrintf(name, 64, "Field_%" PetscInt_FMT, f));
396: PetscCall(PetscStrallocpy(name, &s->fieldNames[f]));
397: PetscCall(PetscSNPrintf(name, 64, "Component_0"));
398: PetscCall(PetscMalloc1(s->numFieldComponents[f], &s->compNames[f]));
399: PetscCall(PetscStrallocpy(name, &s->compNames[f][0]));
400: }
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@
405: PetscSectionGetFieldName - Returns the name of a field in the `PetscSection`
407: Not Collective
409: Input Parameters:
410: + s - the `PetscSection`
411: - field - the field number
413: Output Parameter:
414: . fieldName - the field name
416: Level: intermediate
418: Note:
419: Will error if the field number is out of range
421: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
422: @*/
423: PetscErrorCode PetscSectionGetFieldName(PetscSection s, PetscInt field, const char *fieldName[])
424: {
425: PetscFunctionBegin;
427: PetscAssertPointer(fieldName, 3);
428: PetscSectionCheckValidField(field, s->numFields);
429: *fieldName = s->fieldNames[field];
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: /*@
434: PetscSectionSetFieldName - Sets the name of a field in the `PetscSection`
436: Not Collective
438: Input Parameters:
439: + s - the `PetscSection`
440: . field - the field number
441: - fieldName - the field name
443: Level: intermediate
445: Note:
446: Will error if the field number is out of range
448: .seealso: [PetscSection](ch_petscsection), `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`
449: @*/
450: PetscErrorCode PetscSectionSetFieldName(PetscSection s, PetscInt field, const char fieldName[])
451: {
452: PetscFunctionBegin;
454: if (fieldName) PetscAssertPointer(fieldName, 3);
455: PetscSectionCheckValidField(field, s->numFields);
456: PetscCall(PetscFree(s->fieldNames[field]));
457: PetscCall(PetscStrallocpy(fieldName, &s->fieldNames[field]));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: /*@
462: PetscSectionGetComponentName - Gets the name of a field component in the `PetscSection`
464: Not Collective
466: Input Parameters:
467: + s - the `PetscSection`
468: . field - the field number
469: - comp - the component number
471: Output Parameter:
472: . compName - the component name
474: Level: intermediate
476: Note:
477: Will error if the field or component number do not exist
479: Developer Notes:
480: The function name should have Field in it since they are field components.
482: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
483: `PetscSectionSetComponentName()`, `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
484: @*/
485: PetscErrorCode PetscSectionGetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char *compName[])
486: {
487: PetscFunctionBegin;
489: PetscAssertPointer(compName, 4);
490: PetscSectionCheckValidField(field, s->numFields);
491: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
492: *compName = s->compNames[field][comp];
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: /*@
497: PetscSectionSetComponentName - Sets the name of a field component in the `PetscSection`
499: Not Collective
501: Input Parameters:
502: + s - the `PetscSection`
503: . field - the field number
504: . comp - the component number
505: - compName - the component name
507: Level: advanced
509: Note:
510: Will error if the field or component number do not exist
512: Developer Notes:
513: The function name should have Field in it since they are field components.
515: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetComponentName()`, `PetscSectionSetNumFields()`, `PetscSectionGetNumFields()`,
516: `PetscSectionSetFieldName()`, `PetscSectionGetFieldComponents()`, `PetscSectionSetFieldComponents()`
517: @*/
518: PetscErrorCode PetscSectionSetComponentName(PetscSection s, PetscInt field, PetscInt comp, const char compName[])
519: {
520: PetscFunctionBegin;
522: if (compName) PetscAssertPointer(compName, 4);
523: PetscSectionCheckValidField(field, s->numFields);
524: PetscSectionCheckValidFieldComponent(comp, s->numFieldComponents[field]);
525: PetscCall(PetscFree(s->compNames[field][comp]));
526: PetscCall(PetscStrallocpy(compName, &s->compNames[field][comp]));
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /*@
531: PetscSectionGetFieldComponents - Returns the number of field components for the given field.
533: Not Collective
535: Input Parameters:
536: + s - the `PetscSection`
537: - field - the field number
539: Output Parameter:
540: . numComp - the number of field components
542: Level: advanced
544: Developer Notes:
545: This function is misnamed. There is a Num in `PetscSectionGetNumFields()` but not in this name
547: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldComponents()`, `PetscSectionGetNumFields()`,
548: `PetscSectionSetComponentName()`, `PetscSectionGetComponentName()`
549: @*/
550: PetscErrorCode PetscSectionGetFieldComponents(PetscSection s, PetscInt field, PetscInt *numComp)
551: {
552: PetscFunctionBegin;
554: PetscAssertPointer(numComp, 3);
555: PetscSectionCheckValidField(field, s->numFields);
556: *numComp = s->numFieldComponents[field];
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: /*@
561: PetscSectionSetFieldComponents - Sets the number of field components for the given field.
563: Not Collective
565: Input Parameters:
566: + s - the `PetscSection`
567: . field - the field number
568: - numComp - the number of field components
570: Level: advanced
572: Note:
573: This number can be different than the values set with `PetscSectionSetFieldDof()`. It can be used to indicate the number of
574: components in the field of the underlying physical model which may be different than the number of degrees of freedom needed
575: 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
576: an face based model for velocity (where the velocity normal to the face is stored) there is only 1 dof for each face point.
578: The value set with this function are not needed or used in `PetscSectionSetUp()`.
580: Developer Notes:
581: This function is misnamed. There is a Num in `PetscSectionSetNumFields()` but not in this name
583: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldComponents()`, `PetscSectionSetComponentName()`,
584: `PetscSectionGetComponentName()`, `PetscSectionGetNumFields()`
585: @*/
586: PetscErrorCode PetscSectionSetFieldComponents(PetscSection s, PetscInt field, PetscInt numComp)
587: {
588: PetscFunctionBegin;
590: PetscSectionCheckValidField(field, s->numFields);
591: if (s->compNames) {
592: for (PetscInt c = 0; c < s->numFieldComponents[field]; ++c) PetscCall(PetscFree(s->compNames[field][c]));
593: PetscCall(PetscFree(s->compNames[field]));
594: }
596: s->numFieldComponents[field] = numComp;
597: if (numComp) {
598: PetscCall(PetscMalloc1(numComp, &s->compNames[field]));
599: for (PetscInt c = 0; c < numComp; ++c) {
600: char name[64];
602: PetscCall(PetscSNPrintf(name, 64, "%" PetscInt_FMT, c));
603: PetscCall(PetscStrallocpy(name, &s->compNames[field][c]));
604: }
605: }
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /*@
610: PetscSectionGetChart - Returns the range [`pStart`, `pEnd`) in which points (indices) lie for this `PetscSection` on this MPI process
612: Not Collective
614: Input Parameter:
615: . s - the `PetscSection`
617: Output Parameters:
618: + pStart - the first point
619: - pEnd - one past the last point
621: Level: intermediate
623: Note:
624: The chart may be thought of as the bounds on the points (indices) one may use to index into numerical data that is associated with
625: the `PetscSection` data layout.
627: .seealso: [PetscSection](ch_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 chart may be thought of as the bounds on the points (indices) one may use to index into numerical data that is associated with
652: the `PetscSection` data layout.
654: The charts on different MPI processes may (and often do) overlap
656: If you intend to use `PetscSectionSetNumFields()` it must be called before this call.
658: The chart for all fields created with `PetscSectionSetNumFields()` is the same as this chart.
660: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetChart()`, `PetscSectionCreate()`, `PetscSectionSetNumFields()`
661: @*/
662: PetscErrorCode PetscSectionSetChart(PetscSection s, PetscInt pStart, PetscInt pEnd)
663: {
664: PetscFunctionBegin;
666: PetscCheck(pEnd >= pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Chart pEnd %" PetscInt_FMT " cannot be smaller than chart pStart %" PetscInt_FMT, pEnd, pStart);
667: if (pStart == s->pStart && pEnd == s->pEnd) PetscFunctionReturn(PETSC_SUCCESS);
668: /* Cannot Reset() because it destroys field information */
669: s->setup = PETSC_FALSE;
670: PetscCall(PetscSectionDestroy(&s->bc));
671: PetscCall(PetscFree(s->bcIndices));
672: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
674: s->pStart = pStart;
675: s->pEnd = pEnd;
676: PetscCall(PetscMalloc2(pEnd - pStart, &s->atlasDof, pEnd - pStart, &s->atlasOff));
677: PetscCall(PetscArrayzero(s->atlasDof, pEnd - pStart));
678: for (PetscInt f = 0; f < s->numFields; ++f) PetscCall(PetscSectionSetChart(s->field[f], pStart, pEnd));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: /*@
683: PetscSectionGetPermutation - Returns the permutation of [0, `pEnd` - `pStart`) or `NULL` that was set with `PetscSectionSetPermutation()`
685: Not Collective
687: Input Parameter:
688: . s - the `PetscSection`
690: Output Parameter:
691: . perm - The permutation as an `IS`
693: Level: intermediate
695: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetPermutation()`, `PetscSectionCreate()`
696: @*/
697: PetscErrorCode PetscSectionGetPermutation(PetscSection s, IS *perm)
698: {
699: PetscFunctionBegin;
701: if (perm) {
702: PetscAssertPointer(perm, 2);
703: *perm = s->perm;
704: }
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: /*@
709: PetscSectionSetPermutation - Sets a permutation of the chart for this section, [0, `pEnd` - `pStart`), which determines the order to store the `PetscSection` information
711: Not Collective
713: Input Parameters:
714: + s - the `PetscSection`
715: - perm - the permutation of points
717: Level: intermediate
719: Notes:
720: The permutation must be provided before `PetscSectionSetUp()`.
722: The data in the `PetscSection` are permuted but the access via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` is not changed
724: Compare to `PetscSectionPermute()`
726: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetUp()`, `PetscSectionGetPermutation()`, `PetscSectionPermute()`, `PetscSectionCreate()`
727: @*/
728: PetscErrorCode PetscSectionSetPermutation(PetscSection s, IS perm)
729: {
730: PetscFunctionBegin;
733: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set a permutation after the section is setup");
734: if (s->perm != perm) {
735: PetscCall(ISDestroy(&s->perm));
736: if (perm) {
737: s->perm = perm;
738: PetscCall(PetscObjectReference((PetscObject)s->perm));
739: }
740: }
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: /*@C
745: PetscSectionGetBlockStarts - Returns a table indicating which points start new blocks
747: Not Collective, No Fortran Support
749: Input Parameter:
750: . s - the `PetscSection`
752: Output Parameter:
753: . blockStarts - The `PetscBT` with a 1 for each point that begins a block
755: Notes:
756: The table is on [0, `pEnd` - `pStart`).
758: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
760: Level: intermediate
762: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionSetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
763: @*/
764: PetscErrorCode PetscSectionGetBlockStarts(PetscSection s, PetscBT *blockStarts)
765: {
766: PetscFunctionBegin;
768: if (blockStarts) {
769: PetscAssertPointer(blockStarts, 2);
770: *blockStarts = s->blockStarts;
771: }
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: /*@C
776: PetscSectionSetBlockStarts - Sets a table indicating which points start new blocks
778: Not Collective, No Fortran Support
780: Input Parameters:
781: + s - the `PetscSection`
782: - blockStarts - The `PetscBT` with a 1 for each point that begins a block
784: Level: intermediate
786: Notes:
787: The table is on [0, `pEnd` - `pStart`). PETSc takes ownership of the `PetscBT` when it is passed in and will destroy it. The user should not destroy it.
789: This information is used by `DMCreateMatrix()` to create a variable block size description which is set using `MatSetVariableBlockSizes()`.
791: .seealso: [](sec_scatter), `IS`, `PetscSection`, `PetscSectionGetBlockStarts()`, `PetscSectionCreate()`, `DMCreateMatrix()`, `MatSetVariableBlockSizes()`
792: @*/
793: PetscErrorCode PetscSectionSetBlockStarts(PetscSection s, PetscBT blockStarts)
794: {
795: PetscFunctionBegin;
797: if (s->blockStarts != blockStarts) {
798: PetscCall(PetscBTDestroy(&s->blockStarts));
799: s->blockStarts = blockStarts;
800: }
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: /*@
805: PetscSectionGetPointMajor - Returns the flag for dof ordering, `PETSC_TRUE` if it is point major, `PETSC_FALSE` if it is field major
807: Not Collective
809: Input Parameter:
810: . s - the `PetscSection`
812: Output Parameter:
813: . pm - the flag for point major ordering
815: Level: intermediate
817: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetPointMajor()`
818: @*/
819: PetscErrorCode PetscSectionGetPointMajor(PetscSection s, PetscBool *pm)
820: {
821: PetscFunctionBegin;
823: PetscAssertPointer(pm, 2);
824: *pm = s->pointMajor;
825: PetscFunctionReturn(PETSC_SUCCESS);
826: }
828: /*@
829: PetscSectionSetPointMajor - Sets the flag for dof ordering, `PETSC_TRUE` for point major, otherwise it will be field major
831: Not Collective
833: Input Parameters:
834: + s - the `PetscSection`
835: - pm - the flag for point major ordering
837: Level: intermediate
839: Note:
840: Field-major order is not recommended unless you are managing the entire problem yourself, since many higher-level functions in PETSc depend on point-major order.
842: Point major order means the degrees of freedom are stored as follows
843: .vb
844: all the degrees of freedom for each point are stored contiguously, one point after another (respecting a permutation set with `PetscSectionSetPermutation()`)
845: for each point
846: the degrees of freedom for each field (starting with the unnamed default field) are listed in order by field
847: .ve
849: Field major order means the degrees of freedom are stored as follows
850: .vb
851: all degrees of freedom for each field (including the unnamed default field) are stored contiguously, one field after another
852: for each field (started with unnamed default field)
853: the degrees of freedom for each point are listed in order by point (respecting a permutation set with `PetscSectionSetPermutation()`)
854: .ve
856: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetPointMajor()`, `PetscSectionSetPermutation()`
857: @*/
858: PetscErrorCode PetscSectionSetPointMajor(PetscSection s, PetscBool pm)
859: {
860: PetscFunctionBegin;
862: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the dof ordering after the section is setup");
863: s->pointMajor = pm;
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: PetscSectionGetIncludesConstraints - Returns the flag indicating if constrained dofs were included when computing offsets in the `PetscSection`.
869: The value is set with `PetscSectionSetIncludesConstraints()`
871: Not Collective
873: Input Parameter:
874: . s - the `PetscSection`
876: Output Parameter:
877: . includesConstraints - the flag indicating if constrained dofs were included when computing offsets
879: Level: intermediate
881: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetIncludesConstraints()`
882: @*/
883: PetscErrorCode PetscSectionGetIncludesConstraints(PetscSection s, PetscBool *includesConstraints)
884: {
885: PetscFunctionBegin;
887: PetscAssertPointer(includesConstraints, 2);
888: *includesConstraints = s->includesConstraints;
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: /*@
893: PetscSectionSetIncludesConstraints - Sets the flag indicating if constrained dofs are to be included when computing offsets
895: Not Collective
897: Input Parameters:
898: + s - the `PetscSection`
899: - includesConstraints - the flag indicating if constrained dofs are to be included when computing offsets
901: Level: intermediate
903: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetIncludesConstraints()`
904: @*/
905: PetscErrorCode PetscSectionSetIncludesConstraints(PetscSection s, PetscBool includesConstraints)
906: {
907: PetscFunctionBegin;
909: PetscCheck(!s->setup, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONGSTATE, "Cannot set includesConstraints after the section is set up");
910: s->includesConstraints = includesConstraints;
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: /*@
915: PetscSectionGetDof - Return the total number of degrees of freedom associated with a given point.
917: Not Collective
919: Input Parameters:
920: + s - the `PetscSection`
921: - point - the point
923: Output Parameter:
924: . numDof - the number of dof
926: Level: intermediate
928: Notes:
929: In a global section, this size will be negative for points not owned by this process.
931: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
933: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetDof()`, `PetscSectionCreate()`
934: @*/
935: PetscErrorCode PetscSectionGetDof(PetscSection s, PetscInt point, PetscInt *numDof)
936: {
937: PetscFunctionBeginHot;
939: PetscAssertPointer(numDof, 3);
940: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
941: *numDof = s->atlasDof[point - s->pStart];
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: /*@
946: PetscSectionSetDof - Sets the total number of degrees of freedom associated with a given point.
948: Not Collective
950: Input Parameters:
951: + s - the `PetscSection`
952: . point - the point
953: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
955: Level: intermediate
957: Note:
958: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
960: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionAddDof()`, `PetscSectionCreate()`
961: @*/
962: PetscErrorCode PetscSectionSetDof(PetscSection s, PetscInt point, PetscInt numDof)
963: {
964: PetscFunctionBegin;
966: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
967: s->atlasDof[point - s->pStart] = numDof;
968: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
969: PetscFunctionReturn(PETSC_SUCCESS);
970: }
972: /*@
973: PetscSectionAddDof - Adds to the total number of degrees of freedom associated with a given point.
975: Not Collective
977: Input Parameters:
978: + s - the `PetscSection`
979: . point - the point
980: - numDof - the number of additional dof
982: Level: intermediate
984: Note:
985: This number is for the unnamed default field at the given point plus all degrees of freedom associated with all fields at that point
987: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetDof()`, `PetscSectionCreate()`
988: @*/
989: PetscErrorCode PetscSectionAddDof(PetscSection s, PetscInt point, PetscInt numDof)
990: {
991: PetscFunctionBeginHot;
993: PetscAssert(point >= s->pStart && point < s->pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Section point %" PetscInt_FMT " should be in [%" PetscInt_FMT ", %" PetscInt_FMT ")", point, s->pStart, s->pEnd);
994: PetscCheck(numDof >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numDof %" PetscInt_FMT " should not be negative", numDof);
995: s->atlasDof[point - s->pStart] += numDof;
996: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /*@
1001: PetscSectionGetFieldDof - Return the number of degrees of freedom associated with a field on a given point.
1003: Not Collective
1005: Input Parameters:
1006: + s - the `PetscSection`
1007: . point - the point
1008: - field - the field
1010: Output Parameter:
1011: . numDof - the number of dof
1013: Level: intermediate
1015: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionCreate()`
1016: @*/
1017: PetscErrorCode PetscSectionGetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt *numDof)
1018: {
1019: PetscFunctionBegin;
1021: PetscAssertPointer(numDof, 4);
1022: PetscSectionCheckValidField(field, s->numFields);
1023: PetscCall(PetscSectionGetDof(s->field[field], point, numDof));
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }
1027: /*@
1028: PetscSectionSetFieldDof - Sets the number of degrees of freedom associated with a field on a given point.
1030: Not Collective
1032: Input Parameters:
1033: + s - the `PetscSection`
1034: . point - the point
1035: . field - the field
1036: - numDof - the number of dof, these values may be negative -(dof+1) to indicate they are off process
1038: Level: intermediate
1040: Note:
1041: When setting the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
1042: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1044: This is equivalent to
1045: .vb
1046: PetscSection fs;
1047: PetscSectionGetField(s,field,&fs)
1048: PetscSectionSetDof(fs,numDof)
1049: .ve
1051: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`, `PetscSectionAddDof()`, `PetscSectionSetDof()`
1052: @*/
1053: PetscErrorCode PetscSectionSetFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1054: {
1055: PetscFunctionBegin;
1057: PetscSectionCheckValidField(field, s->numFields);
1058: PetscCall(PetscSectionSetDof(s->field[field], point, numDof));
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: /*@
1063: PetscSectionAddFieldDof - Adds a number of degrees of freedom associated with a field on a given point.
1065: Not Collective
1067: Input Parameters:
1068: + s - the `PetscSection`
1069: . point - the point
1070: . field - the field
1071: - numDof - the number of dof
1073: Level: intermediate
1075: Notes:
1076: When adding to the number of dof for a field at a point one must also ensure the count of the total number of dof at the point (summed over
1077: the fields and the unnamed default field) is correct by also calling `PetscSectionAddDof()` or `PetscSectionSetDof()`
1079: This is equivalent to
1080: .vb
1081: PetscSection fs;
1082: PetscSectionGetField(s,field,&fs)
1083: PetscSectionAddDof(fs,numDof)
1084: .ve
1086: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSetFieldDof()`, `PetscSectionGetFieldDof()`, `PetscSectionCreate()`
1087: @*/
1088: PetscErrorCode PetscSectionAddFieldDof(PetscSection s, PetscInt point, PetscInt field, PetscInt numDof)
1089: {
1090: PetscFunctionBegin;
1092: PetscSectionCheckValidField(field, s->numFields);
1093: PetscCall(PetscSectionAddDof(s->field[field], point, numDof));
1094: PetscFunctionReturn(PETSC_SUCCESS);
1095: }
1097: /*@
1098: PetscSectionGetConstraintDof - Return the number of constrained degrees of freedom associated with a given point.
1100: Not Collective
1102: Input Parameters:
1103: + s - the `PetscSection`
1104: - point - the point
1106: Output Parameter:
1107: . numDof - the number of dof which are fixed by constraints
1109: Level: intermediate
1111: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetDof()`, `PetscSectionSetConstraintDof()`, `PetscSectionCreate()`
1112: @*/
1113: PetscErrorCode PetscSectionGetConstraintDof(PetscSection s, PetscInt point, PetscInt *numDof)
1114: {
1115: PetscFunctionBegin;
1117: PetscAssertPointer(numDof, 3);
1118: if (s->bc) 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](ch_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](ch_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](ch_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](ch_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](ch_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](ch_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](ch_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](ch_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](ch_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](ch_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](ch_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: PetscInt 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)sf)));
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, foff));
1570: PetscCall(PetscSectionGetFieldConstraintDof(gs, p, f, &cdof));
1571: foff = off < 0 ? foff - (dof - cdof) : foff + (dof - cdof);
1572: }
1573: }
1574: for (f = 0; f < numFields; ++f) {
1575: PetscSection gfs = gs->field[f];
1577: PetscCall(PetscSectionSetUpBC(gfs));
1578: 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]));
1579: }
1580: gs->setup = PETSC_TRUE;
1581: PetscCall(PetscSectionViewFromOptions(gs, NULL, "-global_section_view"));
1582: *gsection = gs;
1583: PetscFunctionReturn(PETSC_SUCCESS);
1584: }
1586: /*@
1587: PetscSectionCreateGlobalSectionCensored - Create a `PetscSection` describing the globallayout using
1588: a local (sequential) `PetscSection` on each MPI process and an `PetscSF` describing the section point overlap.
1590: Input Parameters:
1591: + s - The `PetscSection` for the local field layout
1592: . sf - The `PetscSF` describing parallel layout of the section points
1593: . includeConstraints - By default this is `PETSC_FALSE`, meaning that the global vector will not possess constrained dofs
1594: . numExcludes - The number of exclusion ranges, this must have the same value on all MPI processes
1595: - 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
1597: Output Parameter:
1598: . gsection - The `PetscSection` for the global field layout
1600: Level: advanced
1602: Notes:
1603: On each MPI process `gsection` inherits the chart of the `s` on that process.
1605: 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`.
1606: In those locations the value of size is -(size+1) and the value of the offset on the remote process is -(off+1).
1608: This routine augments `PetscSectionCreateGlobalSection()` by allowing one to exclude certain ranges in the chart of the `PetscSection`
1610: Developer Notes:
1611: This is a terrible function name
1613: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`
1614: @*/
1615: PetscErrorCode PetscSectionCreateGlobalSectionCensored(PetscSection s, PetscSF sf, PetscBool includeConstraints, PetscInt numExcludes, const PetscInt excludes[], PetscSection *gsection)
1616: {
1617: const PetscInt *pind = NULL;
1618: PetscInt *neg = NULL, *tmpOff = NULL;
1619: PetscInt pStart, pEnd, p, e, dof, cdof, globalOff = 0, nroots;
1620: PetscInt off;
1622: PetscFunctionBegin;
1625: PetscAssertPointer(gsection, 6);
1626: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), gsection));
1627: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1628: PetscCall(PetscSectionSetChart(*gsection, pStart, pEnd));
1629: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
1630: if (nroots >= 0) {
1631: PetscCheck(nroots >= pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PetscSF nroots %" PetscInt_FMT " < %" PetscInt_FMT " section size", nroots, pEnd - pStart);
1632: PetscCall(PetscCalloc1(nroots, &neg));
1633: if (nroots > pEnd - pStart) {
1634: PetscCall(PetscCalloc1(nroots, &tmpOff));
1635: } else {
1636: tmpOff = &(*gsection)->atlasDof[-pStart];
1637: }
1638: }
1639: /* Mark ghost points with negative dof */
1640: for (p = pStart; p < pEnd; ++p) {
1641: for (e = 0; e < numExcludes; ++e) {
1642: if ((p >= excludes[e * 2 + 0]) && (p < excludes[e * 2 + 1])) {
1643: PetscCall(PetscSectionSetDof(*gsection, p, 0));
1644: break;
1645: }
1646: }
1647: if (e < numExcludes) continue;
1648: PetscCall(PetscSectionGetDof(s, p, &dof));
1649: PetscCall(PetscSectionSetDof(*gsection, p, dof));
1650: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1651: if (!includeConstraints && cdof > 0) PetscCall(PetscSectionSetConstraintDof(*gsection, p, cdof));
1652: if (neg) neg[p] = -(dof + 1);
1653: }
1654: PetscCall(PetscSectionSetUpBC(*gsection));
1655: if (nroots >= 0) {
1656: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1657: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1658: if (nroots > pEnd - pStart) {
1659: for (p = pStart; p < pEnd; ++p) {
1660: if (tmpOff[p] < 0) (*gsection)->atlasDof[p - pStart] = tmpOff[p];
1661: }
1662: }
1663: }
1664: /* Calculate new sizes, get process offset, and calculate point offsets */
1665: if (s->perm) PetscCall(ISGetIndices(s->perm, &pind));
1666: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1667: const PetscInt q = pind ? pind[p] : p;
1669: cdof = (!includeConstraints && s->bc) ? s->bc->atlasDof[q] : 0;
1670: (*gsection)->atlasOff[q] = off;
1671: off += (*gsection)->atlasDof[q] > 0 ? (*gsection)->atlasDof[q] - cdof : 0;
1672: }
1673: PetscCallMPI(MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)s)));
1674: globalOff -= off;
1675: for (p = 0, off = 0; p < pEnd - pStart; ++p) {
1676: (*gsection)->atlasOff[p] += globalOff;
1677: if (neg) neg[p + pStart] = -((*gsection)->atlasOff[p] + 1);
1678: }
1679: if (s->perm) PetscCall(ISRestoreIndices(s->perm, &pind));
1680: /* Put in negative offsets for ghost points */
1681: if (nroots >= 0) {
1682: if (nroots == pEnd - pStart) tmpOff = &(*gsection)->atlasOff[-pStart];
1683: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1684: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, neg, tmpOff, MPI_REPLACE));
1685: if (nroots > pEnd - pStart) {
1686: for (p = pStart; p < pEnd; ++p) {
1687: if (tmpOff[p] < 0) (*gsection)->atlasOff[p - pStart] = tmpOff[p];
1688: }
1689: }
1690: }
1691: if (nroots >= 0 && nroots > pEnd - pStart) PetscCall(PetscFree(tmpOff));
1692: PetscCall(PetscFree(neg));
1693: PetscFunctionReturn(PETSC_SUCCESS);
1694: }
1696: /*@
1697: PetscSectionGetPointLayout - Get a `PetscLayout` for the points with nonzero dof counts of the unnamed default field within this `PetscSection`s local chart
1699: Collective
1701: Input Parameters:
1702: + comm - The `MPI_Comm`
1703: - s - The `PetscSection`
1705: Output Parameter:
1706: . layout - The point layout for the data that defines the section
1708: Level: advanced
1710: Notes:
1711: `PetscSectionGetValueLayout()` provides similar information but counting the total number of degrees of freedom on the MPI process (excluding constrained
1712: degrees of freedom).
1714: This count includes constrained degrees of freedom
1716: This is usually called on the default global section.
1718: Example:
1719: .vb
1720: The chart is [2,5), point 2 has 2 dof, point 3 has 0 dof, point 4 has 1 dof
1721: The local size of the `PetscLayout` is 2 since 2 points have a non-zero number of dof
1722: .ve
1724: Developer Notes:
1725: I find the names of these two functions extremely non-informative
1727: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetValueLayout()`, `PetscSectionCreate()`
1728: @*/
1729: PetscErrorCode PetscSectionGetPointLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1730: {
1731: PetscInt pStart, pEnd, p, localSize = 0;
1733: PetscFunctionBegin;
1734: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1735: for (p = pStart; p < pEnd; ++p) {
1736: PetscInt dof;
1738: PetscCall(PetscSectionGetDof(s, p, &dof));
1739: if (dof >= 0) ++localSize;
1740: }
1741: PetscCall(PetscLayoutCreate(comm, layout));
1742: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1743: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1744: PetscCall(PetscLayoutSetUp(*layout));
1745: PetscFunctionReturn(PETSC_SUCCESS);
1746: }
1748: /*@
1749: PetscSectionGetValueLayout - Get the `PetscLayout` associated with the section dofs of a `PetscSection`
1751: Collective
1753: Input Parameters:
1754: + comm - The `MPI_Comm`
1755: - s - The `PetscSection`
1757: Output Parameter:
1758: . layout - The dof layout for the section
1760: Level: advanced
1762: Notes:
1763: `PetscSectionGetPointLayout()` provides similar information but only counting the number of points with nonzero degrees of freedom and
1764: including the constrained degrees of freedom
1766: This is usually called for the default global section.
1768: Example:
1769: .vb
1770: The chart is [2,5), point 2 has 4 dof (2 constrained), point 3 has 0 dof, point 4 has 1 dof (not constrained)
1771: The local size of the `PetscLayout` is 3 since there are 3 unconstrained degrees of freedom on this MPI process
1772: .ve
1774: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetPointLayout()`, `PetscSectionCreate()`
1775: @*/
1776: PetscErrorCode PetscSectionGetValueLayout(MPI_Comm comm, PetscSection s, PetscLayout *layout)
1777: {
1778: PetscInt pStart, pEnd, p, localSize = 0;
1780: PetscFunctionBegin;
1782: PetscAssertPointer(layout, 3);
1783: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1784: for (p = pStart; p < pEnd; ++p) {
1785: PetscInt dof, cdof;
1787: PetscCall(PetscSectionGetDof(s, p, &dof));
1788: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
1789: if (dof - cdof > 0) localSize += dof - cdof;
1790: }
1791: PetscCall(PetscLayoutCreate(comm, layout));
1792: PetscCall(PetscLayoutSetLocalSize(*layout, localSize));
1793: PetscCall(PetscLayoutSetBlockSize(*layout, 1));
1794: PetscCall(PetscLayoutSetUp(*layout));
1795: PetscFunctionReturn(PETSC_SUCCESS);
1796: }
1798: /*@
1799: PetscSectionGetOffset - Return the offset into an array or `Vec` for the dof associated with the given point.
1801: Not Collective
1803: Input Parameters:
1804: + s - the `PetscSection`
1805: - point - the point
1807: Output Parameter:
1808: . offset - the offset
1810: Level: intermediate
1812: Notes:
1813: In a global section, `offset` will be negative for points not owned by this process.
1815: This is for the unnamed default field in the `PetscSection` not the named fields
1817: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1819: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetPointMajor()`
1820: @*/
1821: PetscErrorCode PetscSectionGetOffset(PetscSection s, PetscInt point, PetscInt *offset)
1822: {
1823: PetscFunctionBegin;
1825: PetscAssertPointer(offset, 3);
1826: 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);
1827: *offset = s->atlasOff[point - s->pStart];
1828: PetscFunctionReturn(PETSC_SUCCESS);
1829: }
1831: /*@
1832: PetscSectionSetOffset - Set the offset into an array or `Vec` for the dof associated with the given point.
1834: Not Collective
1836: Input Parameters:
1837: + s - the `PetscSection`
1838: . point - the point
1839: - offset - the offset, these values may be negative indicating the values are off process
1841: Level: developer
1843: Note:
1844: The user usually does not call this function, but uses `PetscSectionSetUp()`
1846: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1847: @*/
1848: PetscErrorCode PetscSectionSetOffset(PetscSection s, PetscInt point, PetscInt offset)
1849: {
1850: PetscFunctionBegin;
1852: 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);
1853: s->atlasOff[point - s->pStart] = offset;
1854: PetscFunctionReturn(PETSC_SUCCESS);
1855: }
1857: /*@
1858: PetscSectionGetFieldOffset - Return the offset into an array or `Vec` for the field dof associated with the given point.
1860: Not Collective
1862: Input Parameters:
1863: + s - the `PetscSection`
1864: . point - the point
1865: - field - the field
1867: Output Parameter:
1868: . offset - the offset
1870: Level: intermediate
1872: Notes:
1873: In a global section, `offset` will be negative for points not owned by this process.
1875: The `offset` values are different depending on a value set with `PetscSectionSetPointMajor()`
1877: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldPointOffset()`
1878: @*/
1879: PetscErrorCode PetscSectionGetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1880: {
1881: PetscFunctionBegin;
1883: PetscAssertPointer(offset, 4);
1884: PetscSectionCheckValidField(field, s->numFields);
1885: PetscCall(PetscSectionGetOffset(s->field[field], point, offset));
1886: PetscFunctionReturn(PETSC_SUCCESS);
1887: }
1889: /*@
1890: PetscSectionSetFieldOffset - Set the offset into an array or `Vec` for the dof associated with the given field at a point.
1892: Not Collective
1894: Input Parameters:
1895: + s - the `PetscSection`
1896: . point - the point
1897: . field - the field
1898: - offset - the offset, these values may be negative indicating the values are off process
1900: Level: developer
1902: Note:
1903: The user usually does not call this function, but uses `PetscSectionSetUp()`
1905: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetFieldOffset()`, `PetscSectionSetOffset()`, `PetscSectionCreate()`, `PetscSectionSetUp()`
1906: @*/
1907: PetscErrorCode PetscSectionSetFieldOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt offset)
1908: {
1909: PetscFunctionBegin;
1911: PetscSectionCheckValidField(field, s->numFields);
1912: PetscCall(PetscSectionSetOffset(s->field[field], point, offset));
1913: PetscFunctionReturn(PETSC_SUCCESS);
1914: }
1916: /*@
1917: PetscSectionGetFieldPointOffset - Return the offset for the first field dof associated with the given point relative to the offset for that point for the
1918: unnamed default field's first dof
1920: Not Collective
1922: Input Parameters:
1923: + s - the `PetscSection`
1924: . point - the point
1925: - field - the field
1927: Output Parameter:
1928: . offset - the offset
1930: Level: advanced
1932: Note:
1933: This ignores constraints
1935: Example:
1936: .vb
1937: if PetscSectionSetPointMajor(s,PETSC_TRUE)
1938: The unnamed default field has 3 dof at `point`
1939: Field 0 has 2 dof at `point`
1940: Then PetscSectionGetFieldPointOffset(s,point,1,&offset) returns and offset of 5
1941: .ve
1943: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`, `PetscSectionGetFieldOffset()`
1944: @*/
1945: PetscErrorCode PetscSectionGetFieldPointOffset(PetscSection s, PetscInt point, PetscInt field, PetscInt *offset)
1946: {
1947: PetscInt off, foff;
1949: PetscFunctionBegin;
1951: PetscAssertPointer(offset, 4);
1952: PetscSectionCheckValidField(field, s->numFields);
1953: PetscCall(PetscSectionGetOffset(s, point, &off));
1954: PetscCall(PetscSectionGetOffset(s->field[field], point, &foff));
1955: *offset = foff - off;
1956: PetscFunctionReturn(PETSC_SUCCESS);
1957: }
1959: /*@
1960: PetscSectionGetOffsetRange - Return the full range of offsets [`start`, `end`) for a `PetscSection`
1962: Not Collective
1964: Input Parameter:
1965: . s - the `PetscSection`
1967: Output Parameters:
1968: + start - the minimum offset
1969: - end - one more than the maximum offset
1971: Level: intermediate
1973: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetOffset()`, `PetscSectionCreate()`
1974: @*/
1975: PetscErrorCode PetscSectionGetOffsetRange(PetscSection s, PetscInt *start, PetscInt *end)
1976: {
1977: PetscInt os = 0, oe = 0, pStart, pEnd, p;
1979: PetscFunctionBegin;
1981: if (s->atlasOff) {
1982: os = s->atlasOff[0];
1983: oe = s->atlasOff[0];
1984: }
1985: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1986: for (p = 0; p < pEnd - pStart; ++p) {
1987: PetscInt dof = s->atlasDof[p], off = s->atlasOff[p];
1989: if (off >= 0) {
1990: os = PetscMin(os, off);
1991: oe = PetscMax(oe, off + dof);
1992: }
1993: }
1994: if (start) *start = os;
1995: if (end) *end = oe;
1996: PetscFunctionReturn(PETSC_SUCCESS);
1997: }
1999: /*@
2000: PetscSectionCreateSubsection - Create a new, smaller `PetscSection` composed of only selected fields
2002: Collective
2004: Input Parameters:
2005: + s - the `PetscSection`
2006: . len - the number of subfields
2007: - fields - the subfield numbers
2009: Output Parameter:
2010: . subs - the subsection
2012: Level: advanced
2014: Notes:
2015: The chart of `subs` is the same as the chart of `s`
2017: This will error if a fieldnumber is out of range
2019: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2020: @*/
2021: PetscErrorCode PetscSectionCreateSubsection(PetscSection s, PetscInt len, const PetscInt fields[], PetscSection *subs)
2022: {
2023: PetscInt nF, f, c, pStart, pEnd, p, maxCdof = 0;
2025: PetscFunctionBegin;
2026: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2028: PetscAssertPointer(fields, 3);
2029: PetscAssertPointer(subs, 4);
2030: PetscCall(PetscSectionGetNumFields(s, &nF));
2031: 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);
2032: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2033: PetscCall(PetscSectionSetNumFields(*subs, len));
2034: for (f = 0; f < len; ++f) {
2035: const char *name = NULL;
2036: PetscInt numComp = 0;
2037: PetscSectionSym sym;
2039: PetscCall(PetscSectionGetFieldName(s, fields[f], &name));
2040: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2041: PetscCall(PetscSectionGetFieldComponents(s, fields[f], &numComp));
2042: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2043: for (c = 0; c < s->numFieldComponents[fields[f]]; ++c) {
2044: PetscCall(PetscSectionGetComponentName(s, fields[f], c, &name));
2045: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2046: }
2047: PetscCall(PetscSectionGetFieldSym(s, fields[f], &sym));
2048: PetscCall(PetscSectionSetFieldSym(*subs, f, sym));
2049: }
2050: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2051: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2052: for (p = pStart; p < pEnd; ++p) {
2053: PetscInt dof = 0, cdof = 0, fdof = 0, cfdof = 0;
2055: for (f = 0; f < len; ++f) {
2056: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2057: PetscCall(PetscSectionSetFieldDof(*subs, p, f, fdof));
2058: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2059: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, p, f, cfdof));
2060: dof += fdof;
2061: cdof += cfdof;
2062: }
2063: PetscCall(PetscSectionSetDof(*subs, p, dof));
2064: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, p, cdof));
2065: maxCdof = PetscMax(cdof, maxCdof);
2066: }
2067: PetscBT bst, subbst;
2069: PetscCall(PetscSectionGetBlockStarts(s, &bst));
2070: if (bst) {
2071: PetscCall(PetscBTCreate(pEnd - pStart, &subbst));
2072: PetscCall(PetscBTCopy(subbst, pEnd - pStart, bst));
2073: PetscCall(PetscSectionSetBlockStarts(*subs, subbst));
2074: }
2075: PetscCall(PetscSectionSetUp(*subs));
2076: if (maxCdof) {
2077: PetscInt *indices;
2079: PetscCall(PetscMalloc1(maxCdof, &indices));
2080: for (p = pStart; p < pEnd; ++p) {
2081: PetscInt cdof;
2083: PetscCall(PetscSectionGetConstraintDof(*subs, p, &cdof));
2084: if (cdof) {
2085: const PetscInt *oldIndices = NULL;
2086: PetscInt fdof = 0, cfdof = 0, fc, numConst = 0, fOff = 0;
2088: for (f = 0; f < len; ++f) {
2089: PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
2090: PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &cfdof));
2091: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &oldIndices));
2092: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, p, f, oldIndices));
2093: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc] + fOff;
2094: numConst += cfdof;
2095: fOff += fdof;
2096: }
2097: PetscCall(PetscSectionSetConstraintIndices(*subs, p, indices));
2098: }
2099: }
2100: PetscCall(PetscFree(indices));
2101: }
2102: PetscFunctionReturn(PETSC_SUCCESS);
2103: }
2105: /*@
2106: PetscSectionCreateComponentSubsection - Create a new, smaller `PetscSection` composed of only selected components
2108: Collective
2110: Input Parameters:
2111: + s - the `PetscSection`
2112: . len - the number of components
2113: - comps - the component numbers
2115: Output Parameter:
2116: . subs - the subsection
2118: Level: advanced
2120: Notes:
2121: The chart of `subs` is the same as the chart of `s`
2123: This will error if the section has more than one field, or if a component number is out of range
2125: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSupersection()`, `PetscSectionCreate()`
2126: @*/
2127: PetscErrorCode PetscSectionCreateComponentSubsection(PetscSection s, PetscInt len, const PetscInt comps[], PetscSection *subs)
2128: {
2129: PetscSectionSym sym;
2130: const char *name = NULL;
2131: PetscInt Nf, pStart, pEnd;
2133: PetscFunctionBegin;
2134: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2136: PetscAssertPointer(comps, 3);
2137: PetscAssertPointer(subs, 4);
2138: PetscCall(PetscSectionGetNumFields(s, &Nf));
2139: PetscCheck(Nf == 1, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_WRONG, "This method can only handle one field, not %" PetscInt_FMT, Nf);
2140: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2141: PetscCall(PetscSectionSetNumFields(*subs, 1));
2142: PetscCall(PetscSectionGetFieldName(s, 0, &name));
2143: PetscCall(PetscSectionSetFieldName(*subs, 0, name));
2144: PetscCall(PetscSectionSetFieldComponents(*subs, 0, len));
2145: PetscCall(PetscSectionGetFieldSym(s, 0, &sym));
2146: PetscCall(PetscSectionSetFieldSym(*subs, 0, sym));
2147: for (PetscInt c = 0; c < len; ++c) {
2148: PetscCall(PetscSectionGetComponentName(s, 0, comps[c], &name));
2149: PetscCall(PetscSectionSetComponentName(*subs, 0, c, name));
2150: }
2151: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2152: PetscCall(PetscSectionSetChart(*subs, pStart, pEnd));
2153: for (PetscInt p = pStart; p < pEnd; ++p) {
2154: PetscInt dof, cdof, cfdof;
2156: PetscCall(PetscSectionGetDof(s, p, &dof));
2157: if (!dof) continue;
2158: PetscCall(PetscSectionGetFieldConstraintDof(s, p, 0, &cfdof));
2159: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2160: PetscCheck(!cdof && !cfdof, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Component selection does not work with constraints");
2161: PetscCall(PetscSectionSetFieldDof(*subs, p, 0, len));
2162: PetscCall(PetscSectionSetDof(*subs, p, len));
2163: }
2164: PetscCall(PetscSectionSetUp(*subs));
2165: PetscFunctionReturn(PETSC_SUCCESS);
2166: }
2168: /*@
2169: PetscSectionCreateSupersection - Create a new, larger section composed of multiple `PetscSection`s
2171: Collective
2173: Input Parameters:
2174: + s - the input sections
2175: - len - the number of input sections
2177: Output Parameter:
2178: . supers - the supersection
2180: Level: advanced
2182: Notes:
2183: The section offsets now refer to a new, larger vector.
2185: Developer Notes:
2186: Needs to explain how the sections are composed
2188: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubsection()`, `PetscSectionCreate()`
2189: @*/
2190: PetscErrorCode PetscSectionCreateSupersection(PetscSection s[], PetscInt len, PetscSection *supers)
2191: {
2192: PetscInt Nf = 0, f, pStart = PETSC_INT_MAX, pEnd = 0, p, maxCdof = 0, i;
2194: PetscFunctionBegin;
2195: if (!len) PetscFunctionReturn(PETSC_SUCCESS);
2196: for (i = 0; i < len; ++i) {
2197: PetscInt nf, pStarti, pEndi;
2199: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2200: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2201: pStart = PetscMin(pStart, pStarti);
2202: pEnd = PetscMax(pEnd, pEndi);
2203: Nf += nf;
2204: }
2205: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s[0]), supers));
2206: PetscCall(PetscSectionSetNumFields(*supers, Nf));
2207: for (i = 0, f = 0; i < len; ++i) {
2208: PetscInt nf, fi, ci;
2210: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2211: for (fi = 0; fi < nf; ++fi, ++f) {
2212: const char *name = NULL;
2213: PetscInt numComp = 0;
2215: PetscCall(PetscSectionGetFieldName(s[i], fi, &name));
2216: PetscCall(PetscSectionSetFieldName(*supers, f, name));
2217: PetscCall(PetscSectionGetFieldComponents(s[i], fi, &numComp));
2218: PetscCall(PetscSectionSetFieldComponents(*supers, f, numComp));
2219: for (ci = 0; ci < s[i]->numFieldComponents[fi]; ++ci) {
2220: PetscCall(PetscSectionGetComponentName(s[i], fi, ci, &name));
2221: PetscCall(PetscSectionSetComponentName(*supers, f, ci, name));
2222: }
2223: }
2224: }
2225: PetscCall(PetscSectionSetChart(*supers, pStart, pEnd));
2226: for (p = pStart; p < pEnd; ++p) {
2227: PetscInt dof = 0, cdof = 0;
2229: for (i = 0, f = 0; i < len; ++i) {
2230: PetscInt nf, fi, pStarti, pEndi;
2231: PetscInt fdof = 0, cfdof = 0;
2233: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2234: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2235: if ((p < pStarti) || (p >= pEndi)) continue;
2236: for (fi = 0; fi < nf; ++fi, ++f) {
2237: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2238: PetscCall(PetscSectionAddFieldDof(*supers, p, f, fdof));
2239: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2240: if (cfdof) PetscCall(PetscSectionAddFieldConstraintDof(*supers, p, f, cfdof));
2241: dof += fdof;
2242: cdof += cfdof;
2243: }
2244: }
2245: PetscCall(PetscSectionSetDof(*supers, p, dof));
2246: if (cdof) PetscCall(PetscSectionSetConstraintDof(*supers, p, cdof));
2247: maxCdof = PetscMax(cdof, maxCdof);
2248: }
2249: PetscCall(PetscSectionSetUp(*supers));
2250: if (maxCdof) {
2251: PetscInt *indices;
2253: PetscCall(PetscMalloc1(maxCdof, &indices));
2254: for (p = pStart; p < pEnd; ++p) {
2255: PetscInt cdof;
2257: PetscCall(PetscSectionGetConstraintDof(*supers, p, &cdof));
2258: if (cdof) {
2259: PetscInt dof, numConst = 0, fOff = 0;
2261: for (i = 0, f = 0; i < len; ++i) {
2262: const PetscInt *oldIndices = NULL;
2263: PetscInt nf, fi, pStarti, pEndi, fdof, cfdof, fc;
2265: PetscCall(PetscSectionGetNumFields(s[i], &nf));
2266: PetscCall(PetscSectionGetChart(s[i], &pStarti, &pEndi));
2267: if ((p < pStarti) || (p >= pEndi)) continue;
2268: for (fi = 0; fi < nf; ++fi, ++f) {
2269: PetscCall(PetscSectionGetFieldDof(s[i], p, fi, &fdof));
2270: PetscCall(PetscSectionGetFieldConstraintDof(s[i], p, fi, &cfdof));
2271: PetscCall(PetscSectionGetFieldConstraintIndices(s[i], p, fi, &oldIndices));
2272: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] = oldIndices[fc];
2273: PetscCall(PetscSectionSetFieldConstraintIndices(*supers, p, f, &indices[numConst]));
2274: for (fc = 0; fc < cfdof; ++fc) indices[numConst + fc] += fOff;
2275: numConst += cfdof;
2276: }
2277: PetscCall(PetscSectionGetDof(s[i], p, &dof));
2278: fOff += dof;
2279: }
2280: PetscCall(PetscSectionSetConstraintIndices(*supers, p, indices));
2281: }
2282: }
2283: PetscCall(PetscFree(indices));
2284: }
2285: PetscFunctionReturn(PETSC_SUCCESS);
2286: }
2288: static PetscErrorCode PetscSectionCreateSubplexSection_Private(PetscSection s, IS subpointIS, PetscBool renumberPoints, PetscSection *subs)
2289: {
2290: const PetscInt *points = NULL, *indices = NULL;
2291: PetscInt *spoints = NULL, *order = NULL;
2292: PetscInt numFields, f, c, numSubpoints = 0, pStart, pEnd, p, spStart, spEnd, subp;
2294: PetscFunctionBegin;
2297: PetscAssertPointer(subs, 4);
2298: PetscCall(PetscSectionGetNumFields(s, &numFields));
2299: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), subs));
2300: if (numFields) PetscCall(PetscSectionSetNumFields(*subs, numFields));
2301: for (f = 0; f < numFields; ++f) {
2302: const char *name = NULL;
2303: PetscInt numComp = 0;
2305: PetscCall(PetscSectionGetFieldName(s, f, &name));
2306: PetscCall(PetscSectionSetFieldName(*subs, f, name));
2307: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
2308: PetscCall(PetscSectionSetFieldComponents(*subs, f, numComp));
2309: for (c = 0; c < s->numFieldComponents[f]; ++c) {
2310: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
2311: PetscCall(PetscSectionSetComponentName(*subs, f, c, name));
2312: }
2313: }
2314: /* For right now, we do not try to squeeze the subchart */
2315: if (subpointIS) {
2316: PetscCall(ISGetLocalSize(subpointIS, &numSubpoints));
2317: PetscCall(ISGetIndices(subpointIS, &points));
2318: }
2319: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
2320: if (renumberPoints) {
2321: PetscBool sorted;
2323: spStart = 0;
2324: spEnd = numSubpoints;
2325: PetscCall(ISSorted(subpointIS, &sorted));
2326: if (!sorted) {
2327: PetscCall(PetscMalloc2(numSubpoints, &spoints, numSubpoints, &order));
2328: PetscCall(PetscArraycpy(spoints, points, numSubpoints));
2329: for (PetscInt i = 0; i < numSubpoints; ++i) order[i] = i;
2330: PetscCall(PetscSortIntWithArray(numSubpoints, spoints, order));
2331: }
2332: } else {
2333: PetscCall(ISGetMinMax(subpointIS, &spStart, &spEnd));
2334: ++spEnd;
2335: }
2336: PetscCall(PetscSectionSetChart(*subs, spStart, spEnd));
2337: for (p = pStart; p < pEnd; ++p) {
2338: PetscInt dof, cdof, fdof = 0, cfdof = 0;
2340: PetscCall(PetscFindInt(p, numSubpoints, spoints ? spoints : points, &subp));
2341: if (subp < 0) continue;
2342: if (!renumberPoints) subp = p;
2343: else subp = order ? order[subp] : subp;
2344: for (f = 0; f < numFields; ++f) {
2345: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
2346: PetscCall(PetscSectionSetFieldDof(*subs, subp, f, fdof));
2347: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cfdof));
2348: if (cfdof) PetscCall(PetscSectionSetFieldConstraintDof(*subs, subp, f, cfdof));
2349: }
2350: PetscCall(PetscSectionGetDof(s, p, &dof));
2351: PetscCall(PetscSectionSetDof(*subs, subp, dof));
2352: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
2353: if (cdof) PetscCall(PetscSectionSetConstraintDof(*subs, subp, cdof));
2354: }
2355: PetscCall(PetscSectionSetUp(*subs));
2356: /* Change offsets to original offsets */
2357: for (p = pStart; p < pEnd; ++p) {
2358: PetscInt off, foff = 0;
2360: PetscCall(PetscFindInt(p, numSubpoints, spoints ? spoints : points, &subp));
2361: if (subp < 0) continue;
2362: if (!renumberPoints) subp = p;
2363: else subp = order ? order[subp] : subp;
2364: for (f = 0; f < numFields; ++f) {
2365: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
2366: PetscCall(PetscSectionSetFieldOffset(*subs, subp, f, foff));
2367: }
2368: PetscCall(PetscSectionGetOffset(s, p, &off));
2369: PetscCall(PetscSectionSetOffset(*subs, subp, off));
2370: }
2371: /* Copy constraint indices */
2372: for (subp = spStart; subp < spEnd; ++subp) {
2373: PetscInt cdof;
2375: PetscCall(PetscSectionGetConstraintDof(*subs, subp, &cdof));
2376: if (cdof) {
2377: for (f = 0; f < numFields; ++f) {
2378: PetscCall(PetscSectionGetFieldConstraintIndices(s, points[subp - spStart], f, &indices));
2379: PetscCall(PetscSectionSetFieldConstraintIndices(*subs, subp, f, indices));
2380: }
2381: PetscCall(PetscSectionGetConstraintIndices(s, points[subp - spStart], &indices));
2382: PetscCall(PetscSectionSetConstraintIndices(*subs, subp, indices));
2383: }
2384: }
2385: if (subpointIS) PetscCall(ISRestoreIndices(subpointIS, &points));
2386: PetscCall(PetscFree2(spoints, order));
2387: PetscFunctionReturn(PETSC_SUCCESS);
2388: }
2390: /*@
2391: PetscSectionCreateSubmeshSection - Create a new, smaller section with support on the submesh
2393: Collective
2395: Input Parameters:
2396: + s - the `PetscSection`
2397: - subpointIS - a sorted list of points in the original mesh which are in the submesh
2399: Output Parameter:
2400: . subs - the subsection
2402: Level: advanced
2404: Notes:
2405: 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))`
2407: Compare this with `PetscSectionCreateSubdomainSection()` that does not map the points numbers to start at zero but leaves them as before
2409: Developer Notes:
2410: 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`
2412: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubdomainSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2413: @*/
2414: PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection s, IS subpointIS, PetscSection *subs)
2415: {
2416: PetscFunctionBegin;
2417: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointIS, PETSC_TRUE, subs));
2418: PetscFunctionReturn(PETSC_SUCCESS);
2419: }
2421: /*@
2422: PetscSectionCreateSubdomainSection - Create a new, smaller section with support on a subdomain of the mesh
2424: Collective
2426: Input Parameters:
2427: + s - the `PetscSection`
2428: - subpointMap - a sorted list of points in the original mesh which are in the subdomain
2430: Output Parameter:
2431: . subs - the subsection
2433: Level: advanced
2435: Notes:
2436: 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`
2437: is `[min(subpointMap),max(subpointMap)+1)`
2439: Compare this with `PetscSectionCreateSubmeshSection()` that maps the point numbers to start at zero
2441: Developer Notes:
2442: 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`
2444: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreateSubmeshSection()`, `PetscSectionCreateSubsection()`, `DMPlexGetSubpointMap()`, `PetscSectionCreate()`
2445: @*/
2446: PetscErrorCode PetscSectionCreateSubdomainSection(PetscSection s, IS subpointMap, PetscSection *subs)
2447: {
2448: PetscFunctionBegin;
2449: PetscCall(PetscSectionCreateSubplexSection_Private(s, subpointMap, PETSC_FALSE, subs));
2450: PetscFunctionReturn(PETSC_SUCCESS);
2451: }
2453: static PetscErrorCode PetscSectionView_ASCII(PetscSection s, PetscViewer viewer)
2454: {
2455: PetscMPIInt rank;
2457: PetscFunctionBegin;
2458: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2459: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2460: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2461: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
2462: if (s->bc && s->bc->atlasDof[p] > 0) {
2463: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT " constrained", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2464: if (s->bcIndices) {
2465: for (PetscInt b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2466: }
2467: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2468: } else {
2469: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT "\n", p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2470: }
2471: }
2472: PetscCall(PetscViewerFlush(viewer));
2473: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2474: if (s->sym) {
2475: PetscCall(PetscViewerASCIIPushTab(viewer));
2476: PetscCall(PetscSectionSymView(s->sym, viewer));
2477: PetscCall(PetscViewerASCIIPopTab(viewer));
2478: }
2479: PetscFunctionReturn(PETSC_SUCCESS);
2480: }
2482: /*@
2483: PetscSectionViewFromOptions - View the `PetscSection` based on values in the options database
2485: Collective
2487: Input Parameters:
2488: + A - the `PetscSection` object to view
2489: . obj - Optional object that provides the options prefix used for the options
2490: - name - command line option
2492: Options Database Key:
2493: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
2495: Level: intermediate
2497: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionView`, `PetscObjectViewFromOptions()`, `PetscSectionCreate()`, `PetscSectionView()`
2498: @*/
2499: PetscErrorCode PetscSectionViewFromOptions(PetscSection A, PetscObject obj, const char name[])
2500: {
2501: PetscFunctionBegin;
2503: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
2504: PetscFunctionReturn(PETSC_SUCCESS);
2505: }
2507: /*@
2508: PetscSectionView - Views a `PetscSection`
2510: Collective
2512: Input Parameters:
2513: + s - the `PetscSection` object to view
2514: - viewer - the viewer
2516: Level: beginner
2518: Note:
2519: `PetscSectionView()`, when viewer is of type `PETSCVIEWERHDF5`, only saves
2520: distribution independent data, such as dofs, offsets, constraint dofs,
2521: and constraint indices. Points that have negative dofs, for instance,
2522: are not saved as they represent points owned by other processes.
2523: Point numbering and rank assignment is currently not stored.
2524: The saved section can be loaded with `PetscSectionLoad()`.
2526: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionLoad()`, `PetscViewer`
2527: @*/
2528: PetscErrorCode PetscSectionView(PetscSection s, PetscViewer viewer)
2529: {
2530: PetscBool isascii, ishdf5;
2531: PetscInt f;
2533: PetscFunctionBegin;
2535: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2537: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2538: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2539: if (isascii) {
2540: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)s, viewer));
2541: if (s->numFields) {
2542: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " fields\n", s->numFields));
2543: for (f = 0; f < s->numFields; ++f) {
2544: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " \"%s\" with %" PetscInt_FMT " components\n", f, s->fieldNames[f], s->numFieldComponents[f]));
2545: PetscCall(PetscSectionView_ASCII(s->field[f], viewer));
2546: }
2547: } else {
2548: PetscCall(PetscSectionView_ASCII(s, viewer));
2549: }
2550: } else if (ishdf5) {
2551: #if PetscDefined(HAVE_HDF5)
2552: PetscCall(PetscSectionView_HDF5_Internal(s, viewer));
2553: #else
2554: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2555: #endif
2556: }
2557: PetscFunctionReturn(PETSC_SUCCESS);
2558: }
2560: /*@
2561: PetscSectionLoad - Loads a `PetscSection`
2563: Collective
2565: Input Parameters:
2566: + s - the `PetscSection` object to load
2567: - viewer - the viewer
2569: Level: beginner
2571: Note:
2572: `PetscSectionLoad()`, when viewer is of type `PETSCVIEWERHDF5`, loads
2573: a section saved with `PetscSectionView()`. The number of processes
2574: used here (N) does not need to be the same as that used when saving.
2575: After calling this function, the chart of s on rank i will be set
2576: to [0, E_i), where \sum_{i=0}^{N-1}E_i equals to the total number of
2577: saved section points.
2579: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionDestroy()`, `PetscSectionView()`
2580: @*/
2581: PetscErrorCode PetscSectionLoad(PetscSection s, PetscViewer viewer)
2582: {
2583: PetscBool ishdf5;
2585: PetscFunctionBegin;
2588: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2589: PetscCheck(ishdf5, PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "Viewer type %s not yet supported for PetscSection loading", ((PetscObject)viewer)->type_name);
2590: #if PetscDefined(HAVE_HDF5)
2591: PetscCall(PetscSectionLoad_HDF5_Internal(s, viewer));
2592: PetscFunctionReturn(PETSC_SUCCESS);
2593: #else
2594: SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
2595: #endif
2596: }
2598: static inline PetscErrorCode PrintArrayElement(void *array, PetscDataType data_type, PetscCount index, PetscViewer viewer)
2599: {
2600: PetscFunctionBeginUser;
2601: switch (data_type) {
2602: case PETSC_INT: {
2603: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt_FMT, ((PetscInt *)array)[index]));
2604: break;
2605: }
2606: case PETSC_INT32: {
2607: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt32_FMT, ((PetscInt32 *)array)[index]));
2608: break;
2609: }
2610: case PETSC_INT64: {
2611: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscInt64_FMT, ((PetscInt64 *)array)[index]));
2612: break;
2613: }
2614: case PETSC_COUNT: {
2615: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %2" PetscCount_FMT, ((PetscCount *)array)[index]));
2616: break;
2617: }
2618: // PETSC_SCALAR is set to the appropriate type
2619: case PETSC_DOUBLE: {
2620: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", ((double *)array)[index]));
2621: break;
2622: }
2623: case PETSC_FLOAT: {
2624: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((float *)array)[index]));
2625: break;
2626: }
2627: #if defined(PETSC_USE_REAL___FLOAT128)
2628: case PETSC___FLOAT128: {
2629: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((PetscReal *)array)[index]));
2630: break;
2631: }
2632: #endif
2633: #if defined(PETSC_USE_REAL___FP16)
2634: case PETSC___FP16: {
2635: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)((PetscReal *)array)[index]));
2636: break;
2637: }
2638: #endif
2639: #if defined(PETSC_HAVE_COMPLEX)
2640: case PETSC_COMPLEX: {
2641: PetscComplex v = ((PetscComplex *)array)[index];
2642: if (PetscImaginaryPartComplex(v) > 0.0) {
2643: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g + %g i", (double)PetscRealPartComplex(v), (double)PetscImaginaryPartComplex(v)));
2644: } else if (PetscImaginaryPartComplex(v) < 0.0) {
2645: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g - %g i", (double)PetscRealPartComplex(v), (double)(-PetscImaginaryPartComplex(v))));
2646: } else {
2647: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %g", (double)PetscRealPartComplex(v)));
2648: }
2649: break;
2650: }
2651: #endif
2652: default:
2653: SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "PetscDataType %d (%s) not supported", data_type, PetscDataTypes[data_type]);
2654: }
2655: PetscFunctionReturn(PETSC_SUCCESS);
2656: }
2658: PetscErrorCode PetscSectionArrayView_ASCII_Internal(PetscSection s, void *array, PetscDataType data_type, PetscViewer viewer)
2659: {
2660: PetscInt i;
2661: PetscMPIInt rank;
2663: PetscFunctionBegin;
2664: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
2665: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2666: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Process %d:\n", rank));
2667: for (PetscInt p = 0; p < s->pEnd - s->pStart; ++p) {
2668: if (s->bc && (s->bc->atlasDof[p] > 0)) {
2669: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2670: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) PetscCall(PrintArrayElement(array, data_type, i, viewer));
2671: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " constrained"));
2672: for (PetscInt b = 0; b < s->bc->atlasDof[p]; ++b) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, s->bcIndices[s->bc->atlasOff[p] + b]));
2673: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2674: } else {
2675: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%4" PetscInt_FMT ") dof %2" PetscInt_FMT " offset %3" PetscInt_FMT, p + s->pStart, s->atlasDof[p], s->atlasOff[p]));
2676: for (i = s->atlasOff[p]; i < s->atlasOff[p] + s->atlasDof[p]; ++i) PetscCall(PrintArrayElement(array, data_type, i, viewer));
2677: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2678: }
2679: }
2680: PetscCall(PetscViewerFlush(viewer));
2681: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2682: PetscFunctionReturn(PETSC_SUCCESS);
2683: }
2685: /*@
2686: PetscSectionArrayView - View an array, using the section to structure the values
2688: Collective
2690: Input Parameters:
2691: + s - the organizing `PetscSection`
2692: . array - the array of values
2693: . data_type - the `PetscDataType` of the array
2694: - viewer - the `PetscViewer`
2696: Level: developer
2698: .seealso: `PetscSection`, `PetscViewer`, `PetscSectionCreate()`, `VecSetValuesSection()`, `PetscSectionVecView()`
2699: @*/
2700: PetscErrorCode PetscSectionArrayView(PetscSection s, void *array, PetscDataType data_type, PetscViewer viewer)
2701: {
2702: PetscBool isascii;
2703: PetscInt f;
2705: PetscFunctionBegin;
2707: if (!array) {
2708: PetscInt size;
2709: PetscCall(PetscSectionGetStorageSize(s, &size));
2710: PetscCheck(size == 0, PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_SIZ, "NULL array passed, but section's storage size is non-zero");
2711: } else PetscAssertPointer(array, 2);
2712: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)s), &viewer));
2714: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2715: if (isascii) {
2716: if (s->numFields) {
2717: PetscCall(PetscViewerASCIIPrintf(viewer, "Array with %" PetscInt_FMT " fields\n", s->numFields));
2718: for (f = 0; f < s->numFields; ++f) {
2719: PetscCall(PetscViewerASCIIPrintf(viewer, " field %" PetscInt_FMT " with %" PetscInt_FMT " components\n", f, s->numFieldComponents[f]));
2720: PetscCall(PetscSectionArrayView_ASCII_Internal(s->field[f], array, data_type, viewer));
2721: }
2722: } else {
2723: PetscCall(PetscSectionArrayView_ASCII_Internal(s, array, data_type, viewer));
2724: }
2725: }
2726: PetscFunctionReturn(PETSC_SUCCESS);
2727: }
2729: /*@
2730: PetscSectionResetClosurePermutation - Remove any existing closure permutation
2732: Input Parameter:
2733: . section - The `PetscSection`
2735: Level: intermediate
2737: .seealso: `PetscSectionSetClosurePermutation()`, `PetscSectionSetClosureIndex()`, `PetscSectionReset()`
2738: @*/
2739: PetscErrorCode PetscSectionResetClosurePermutation(PetscSection section)
2740: {
2741: PetscSectionClosurePermVal clVal;
2743: PetscFunctionBegin;
2744: if (!section->clHash) PetscFunctionReturn(PETSC_SUCCESS);
2745: kh_foreach_value(section->clHash, clVal, {
2746: PetscCall(PetscFree(clVal.perm));
2747: PetscCall(PetscFree(clVal.invPerm));
2748: });
2749: kh_destroy(ClPerm, section->clHash);
2750: section->clHash = NULL;
2751: PetscFunctionReturn(PETSC_SUCCESS);
2752: }
2754: /*@
2755: PetscSectionReset - Frees all section data, the section is then as if `PetscSectionCreate()` had just been called.
2757: Not Collective
2759: Input Parameter:
2760: . s - the `PetscSection`
2762: Level: beginner
2764: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`
2765: @*/
2766: PetscErrorCode PetscSectionReset(PetscSection s)
2767: {
2768: PetscInt f, c;
2770: PetscFunctionBegin;
2772: for (f = 0; f < s->numFields; ++f) {
2773: PetscCall(PetscSectionDestroy(&s->field[f]));
2774: PetscCall(PetscFree(s->fieldNames[f]));
2775: for (c = 0; c < s->numFieldComponents[f]; ++c) PetscCall(PetscFree(s->compNames[f][c]));
2776: PetscCall(PetscFree(s->compNames[f]));
2777: }
2778: PetscCall(PetscFree(s->numFieldComponents));
2779: PetscCall(PetscFree(s->fieldNames));
2780: PetscCall(PetscFree(s->compNames));
2781: PetscCall(PetscFree(s->field));
2782: PetscCall(PetscSectionDestroy(&s->bc));
2783: PetscCall(PetscFree(s->bcIndices));
2784: PetscCall(PetscFree2(s->atlasDof, s->atlasOff));
2785: PetscCall(PetscSectionDestroy(&s->clSection));
2786: PetscCall(ISDestroy(&s->clPoints));
2787: PetscCall(ISDestroy(&s->perm));
2788: PetscCall(PetscBTDestroy(&s->blockStarts));
2789: PetscCall(PetscSectionResetClosurePermutation(s));
2790: PetscCall(PetscSectionSymDestroy(&s->sym));
2791: PetscCall(PetscSectionDestroy(&s->clSection));
2792: PetscCall(ISDestroy(&s->clPoints));
2793: PetscCall(PetscSectionInvalidateMaxDof_Internal(s));
2794: s->pStart = -1;
2795: s->pEnd = -1;
2796: s->maxDof = 0;
2797: s->setup = PETSC_FALSE;
2798: s->numFields = 0;
2799: s->clObj = NULL;
2800: PetscFunctionReturn(PETSC_SUCCESS);
2801: }
2803: /*@
2804: PetscSectionDestroy - Frees a `PetscSection`
2806: Not Collective
2808: Input Parameter:
2809: . s - the `PetscSection`
2811: Level: beginner
2813: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionCreate()`, `PetscSectionReset()`
2814: @*/
2815: PetscErrorCode PetscSectionDestroy(PetscSection *s)
2816: {
2817: PetscFunctionBegin;
2818: if (!*s) PetscFunctionReturn(PETSC_SUCCESS);
2820: if (--((PetscObject)*s)->refct > 0) {
2821: *s = NULL;
2822: PetscFunctionReturn(PETSC_SUCCESS);
2823: }
2824: PetscCall(PetscSectionReset(*s));
2825: PetscCall(PetscHeaderDestroy(s));
2826: PetscFunctionReturn(PETSC_SUCCESS);
2827: }
2829: static PetscErrorCode VecIntGetValuesSection_Private(const PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt **values)
2830: {
2831: const PetscInt p = point - s->pStart;
2833: PetscFunctionBegin;
2835: *values = &baseArray[s->atlasOff[p]];
2836: PetscFunctionReturn(PETSC_SUCCESS);
2837: }
2839: static PetscErrorCode VecIntSetValuesSection_Private(PetscInt *baseArray, PetscSection s, PetscInt point, const PetscInt values[], InsertMode mode)
2840: {
2841: PetscInt *array;
2842: const PetscInt p = point - s->pStart;
2843: const PetscInt orientation = 0; /* Needs to be included for use in closure operations */
2844: PetscInt cDim = 0;
2846: PetscFunctionBegin;
2848: PetscCall(PetscSectionGetConstraintDof(s, p, &cDim));
2849: array = &baseArray[s->atlasOff[p]];
2850: if (!cDim) {
2851: if (orientation >= 0) {
2852: const PetscInt dim = s->atlasDof[p];
2853: PetscInt i;
2855: if (mode == INSERT_VALUES) {
2856: for (i = 0; i < dim; ++i) array[i] = values ? values[i] : i;
2857: } else {
2858: for (i = 0; i < dim; ++i) array[i] += values[i];
2859: }
2860: } else {
2861: PetscInt offset = 0;
2862: PetscInt j = -1, field, i;
2864: for (field = 0; field < s->numFields; ++field) {
2865: const PetscInt dim = s->field[field]->atlasDof[p];
2867: for (i = dim - 1; i >= 0; --i) array[++j] = values ? values[i + offset] : i + offset;
2868: offset += dim;
2869: }
2870: }
2871: } else {
2872: if (orientation >= 0) {
2873: const PetscInt dim = s->atlasDof[p];
2874: PetscInt cInd = 0, i;
2875: const PetscInt *cDof;
2877: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2878: if (mode == INSERT_VALUES) {
2879: for (i = 0; i < dim; ++i) {
2880: if ((cInd < cDim) && (i == cDof[cInd])) {
2881: ++cInd;
2882: continue;
2883: }
2884: array[i] = values ? values[i] : i;
2885: }
2886: } else {
2887: for (i = 0; i < dim; ++i) {
2888: if ((cInd < cDim) && (i == cDof[cInd])) {
2889: ++cInd;
2890: continue;
2891: }
2892: array[i] += values[i];
2893: }
2894: }
2895: } else {
2896: const PetscInt *cDof;
2897: PetscInt offset = 0;
2898: PetscInt cOffset = 0;
2899: PetscInt j = 0, field;
2901: PetscCall(PetscSectionGetConstraintIndices(s, point, &cDof));
2902: for (field = 0; field < s->numFields; ++field) {
2903: const PetscInt dim = s->field[field]->atlasDof[p]; /* PetscSectionGetFieldDof() */
2904: const PetscInt tDim = s->field[field]->bc->atlasDof[p]; /* PetscSectionGetFieldConstraintDof() */
2905: const PetscInt sDim = dim - tDim;
2906: PetscInt cInd = 0, i, k;
2908: for (i = 0, k = dim + offset - 1; i < dim; ++i, ++j, --k) {
2909: if ((cInd < sDim) && (j == cDof[cInd + cOffset])) {
2910: ++cInd;
2911: continue;
2912: }
2913: array[j] = values ? values[k] : k;
2914: }
2915: offset += dim;
2916: cOffset += dim - tDim;
2917: }
2918: }
2919: }
2920: PetscFunctionReturn(PETSC_SUCCESS);
2921: }
2923: /*@
2924: PetscSectionHasConstraints - Determine whether a `PetscSection` has constrained dofs
2926: Not Collective
2928: Input Parameter:
2929: . s - The `PetscSection`
2931: Output Parameter:
2932: . hasConstraints - flag indicating that the section has constrained dofs
2934: Level: intermediate
2936: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2937: @*/
2938: PetscErrorCode PetscSectionHasConstraints(PetscSection s, PetscBool *hasConstraints)
2939: {
2940: PetscFunctionBegin;
2942: PetscAssertPointer(hasConstraints, 2);
2943: *hasConstraints = s->bc ? PETSC_TRUE : PETSC_FALSE;
2944: PetscFunctionReturn(PETSC_SUCCESS);
2945: }
2947: /*@C
2948: PetscSectionGetConstraintIndices - Get the point dof numbers, in [0, dof), which are constrained for a given point
2950: Not Collective
2952: Input Parameters:
2953: + s - The `PetscSection`
2954: - point - The point
2956: Output Parameter:
2957: . indices - The constrained dofs
2959: Level: intermediate
2961: Fortran Notes:
2962: Use `PetscSectionRestoreConstraintIndices()` when the indices are no longer needed
2964: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2965: @*/
2966: PetscErrorCode PetscSectionGetConstraintIndices(PetscSection s, PetscInt point, const PetscInt *indices[])
2967: {
2968: PetscFunctionBegin;
2970: if (s->bc) PetscCall(VecIntGetValuesSection_Private(s->bcIndices, s->bc, point, indices));
2971: else *indices = NULL;
2972: PetscFunctionReturn(PETSC_SUCCESS);
2973: }
2975: /*@
2976: PetscSectionSetConstraintIndices - Set the point dof numbers, in [0, dof), which are constrained
2978: Not Collective
2980: Input Parameters:
2981: + s - The `PetscSection`
2982: . point - The point
2983: - indices - The constrained dofs
2985: Level: intermediate
2987: .seealso: [PetscSection](ch_petscsection), `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
2988: @*/
2989: PetscErrorCode PetscSectionSetConstraintIndices(PetscSection s, PetscInt point, const PetscInt indices[])
2990: {
2991: PetscFunctionBegin;
2993: if (s->bc) {
2994: const PetscInt dof = s->atlasDof[point];
2995: const PetscInt cdof = s->bc->atlasDof[point];
2996: if (indices)
2997: for (PetscInt d = 0; d < cdof; ++d)
2998: 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]);
2999: PetscCall(VecIntSetValuesSection_Private(s->bcIndices, s->bc, point, indices, INSERT_VALUES));
3000: }
3001: PetscFunctionReturn(PETSC_SUCCESS);
3002: }
3004: /*@C
3005: PetscSectionGetFieldConstraintIndices - Get the field dof numbers, in [0, fdof), which are constrained
3007: Not Collective
3009: Input Parameters:
3010: + s - The `PetscSection`
3011: . field - The field number
3012: - point - The point
3014: Output Parameter:
3015: . indices - The constrained dofs sorted in ascending order, the length is returned by `PetscSectionGetConstraintDof()`.
3017: Level: intermediate
3019: Fortran Notes:
3020: Use `PetscSectionRestoreFieldConstraintIndices()` to restore the indices when no longer needed
3022: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetFieldConstraintIndices()`, `PetscSectionGetConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
3023: @*/
3024: PetscErrorCode PetscSectionGetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt *indices[])
3025: {
3026: PetscFunctionBegin;
3028: PetscAssertPointer(indices, 4);
3029: PetscSectionCheckValidField(field, s->numFields);
3030: PetscCall(PetscSectionGetConstraintIndices(s->field[field], point, indices));
3031: PetscFunctionReturn(PETSC_SUCCESS);
3032: }
3034: /*@
3035: PetscSectionSetFieldConstraintIndices - Set the field dof numbers, in [0, fdof), which are constrained
3037: Not Collective
3039: Input Parameters:
3040: + s - The `PetscSection`
3041: . point - The point
3042: . field - The field number
3043: - indices - The constrained dofs
3045: Level: intermediate
3047: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetConstraintIndices()`, `PetscSectionGetFieldConstraintIndices()`, `PetscSectionGetConstraintDof()`, `PetscSection`
3048: @*/
3049: PetscErrorCode PetscSectionSetFieldConstraintIndices(PetscSection s, PetscInt point, PetscInt field, const PetscInt indices[])
3050: {
3051: PetscFunctionBegin;
3053: PetscSectionCheckValidField(field, s->numFields);
3054: PetscCall(PetscSectionSetConstraintIndices(s->field[field], point, indices));
3055: PetscFunctionReturn(PETSC_SUCCESS);
3056: }
3058: /*@
3059: PetscSectionPermute - Reorder the section according to the input point permutation
3061: Collective
3063: Input Parameters:
3064: + section - The `PetscSection` object
3065: - permutation - The point permutation, old point p becomes new point perm[p]
3067: Output Parameter:
3068: . sectionNew - The permuted `PetscSection`
3070: Level: intermediate
3072: Note:
3073: The data and the access to the data via `PetscSectionGetFieldOffset()` and `PetscSectionGetOffset()` are both changed in `sectionNew`
3075: Compare to `PetscSectionSetPermutation()`
3077: .seealso: [PetscSection](ch_petscsection), `IS`, `PetscSection`, `MatPermute()`, `PetscSectionSetPermutation()`
3078: @*/
3079: PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
3080: {
3081: PetscSection s = section, sNew;
3082: const PetscInt *perm;
3083: PetscInt numFields, f, c, numPoints, pStart, pEnd, p;
3085: PetscFunctionBegin;
3088: PetscAssertPointer(sectionNew, 3);
3089: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)s), &sNew));
3090: PetscCall(PetscSectionGetNumFields(s, &numFields));
3091: if (numFields) PetscCall(PetscSectionSetNumFields(sNew, numFields));
3092: for (f = 0; f < numFields; ++f) {
3093: const char *name;
3094: PetscInt numComp;
3096: PetscCall(PetscSectionGetFieldName(s, f, &name));
3097: PetscCall(PetscSectionSetFieldName(sNew, f, name));
3098: PetscCall(PetscSectionGetFieldComponents(s, f, &numComp));
3099: PetscCall(PetscSectionSetFieldComponents(sNew, f, numComp));
3100: for (c = 0; c < s->numFieldComponents[f]; ++c) {
3101: PetscCall(PetscSectionGetComponentName(s, f, c, &name));
3102: PetscCall(PetscSectionSetComponentName(sNew, f, c, name));
3103: }
3104: }
3105: PetscCall(ISGetLocalSize(permutation, &numPoints));
3106: PetscCall(ISGetIndices(permutation, &perm));
3107: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
3108: PetscCall(PetscSectionSetChart(sNew, pStart, pEnd));
3109: PetscCheck(numPoints >= pEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %" PetscInt_FMT " is less than largest Section point %" PetscInt_FMT, numPoints, pEnd);
3110: for (p = pStart; p < pEnd; ++p) {
3111: PetscInt dof, cdof;
3113: PetscCall(PetscSectionGetDof(s, p, &dof));
3114: PetscCall(PetscSectionSetDof(sNew, perm[p], dof));
3115: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
3116: if (cdof) PetscCall(PetscSectionSetConstraintDof(sNew, perm[p], cdof));
3117: for (f = 0; f < numFields; ++f) {
3118: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
3119: PetscCall(PetscSectionSetFieldDof(sNew, perm[p], f, dof));
3120: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
3121: if (cdof) PetscCall(PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof));
3122: }
3123: }
3124: PetscCall(PetscSectionSetUp(sNew));
3125: for (p = pStart; p < pEnd; ++p) {
3126: const PetscInt *cind;
3127: PetscInt cdof;
3129: PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
3130: if (cdof) {
3131: PetscCall(PetscSectionGetConstraintIndices(s, p, &cind));
3132: PetscCall(PetscSectionSetConstraintIndices(sNew, perm[p], cind));
3133: }
3134: for (f = 0; f < numFields; ++f) {
3135: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
3136: if (cdof) {
3137: PetscCall(PetscSectionGetFieldConstraintIndices(s, p, f, &cind));
3138: PetscCall(PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind));
3139: }
3140: }
3141: }
3142: PetscCall(ISRestoreIndices(permutation, &perm));
3143: *sectionNew = sNew;
3144: PetscFunctionReturn(PETSC_SUCCESS);
3145: }
3147: /*@
3148: PetscSectionSetClosureIndex - Create an internal data structure to speed up closure queries.
3150: Collective
3152: Input Parameters:
3153: + section - The `PetscSection`
3154: . obj - A `PetscObject` which serves as the key for this index
3155: . clSection - `PetscSection` giving the size of the closure of each point
3156: - clPoints - `IS` giving the points in each closure
3158: Level: advanced
3160: Note:
3161: This function creates an internal map from each point to its closure. We compress out closure points with no dofs in this section.
3163: Developer Notes:
3164: The information provided here is completely opaque
3166: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`
3167: @*/
3168: PetscErrorCode PetscSectionSetClosureIndex(PetscSection section, PetscObject obj, PetscSection clSection, IS clPoints)
3169: {
3170: PetscFunctionBegin;
3174: if (section->clObj != obj) PetscCall(PetscSectionResetClosurePermutation(section));
3175: section->clObj = obj;
3176: PetscCall(PetscObjectReference((PetscObject)clSection));
3177: PetscCall(PetscObjectReference((PetscObject)clPoints));
3178: PetscCall(PetscSectionDestroy(§ion->clSection));
3179: PetscCall(ISDestroy(§ion->clPoints));
3180: section->clSection = clSection;
3181: section->clPoints = clPoints;
3182: PetscFunctionReturn(PETSC_SUCCESS);
3183: }
3185: /*@
3186: PetscSectionGetClosureIndex - Get the cache of points in the closure of each point in the section set with `PetscSectionSetClosureIndex()`
3188: Collective
3190: Input Parameters:
3191: + section - The `PetscSection`
3192: - obj - A `PetscObject` which serves as the key for this index
3194: Output Parameters:
3195: + clSection - `PetscSection` giving the size of the closure of each point
3196: - clPoints - `IS` giving the points in each closure
3198: Level: advanced
3200: .seealso: [PetscSection](ch_petscsection), `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3201: @*/
3202: PetscErrorCode PetscSectionGetClosureIndex(PetscSection section, PetscObject obj, PetscSection *clSection, IS *clPoints)
3203: {
3204: PetscFunctionBegin;
3205: if (section->clObj == obj) {
3206: if (clSection) *clSection = section->clSection;
3207: if (clPoints) *clPoints = section->clPoints;
3208: } else {
3209: if (clSection) *clSection = NULL;
3210: if (clPoints) *clPoints = NULL;
3211: }
3212: PetscFunctionReturn(PETSC_SUCCESS);
3213: }
3215: PetscErrorCode PetscSectionSetClosurePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, PetscCopyMode mode, PetscInt *clPerm)
3216: {
3217: khiter_t iter;
3218: int new_entry;
3219: PetscSectionClosurePermKey key = {depth, clSize};
3220: PetscSectionClosurePermVal *val;
3222: PetscFunctionBegin;
3223: if (section->clObj != obj) {
3224: PetscCall(PetscSectionDestroy(§ion->clSection));
3225: PetscCall(ISDestroy(§ion->clPoints));
3226: }
3227: section->clObj = obj;
3228: if (!section->clHash) PetscCall(PetscClPermCreate(§ion->clHash));
3229: iter = kh_put(ClPerm, section->clHash, key, &new_entry);
3230: val = &kh_val(section->clHash, iter);
3231: if (!new_entry) {
3232: PetscCall(PetscFree(val->perm));
3233: PetscCall(PetscFree(val->invPerm));
3234: }
3235: if (mode == PETSC_COPY_VALUES) {
3236: PetscCall(PetscMalloc1(clSize, &val->perm));
3237: PetscCall(PetscArraycpy(val->perm, clPerm, clSize));
3238: } else if (mode == PETSC_OWN_POINTER) {
3239: val->perm = clPerm;
3240: } else SETERRQ(PetscObjectComm(obj), PETSC_ERR_SUP, "Do not support borrowed arrays");
3241: PetscCall(PetscMalloc1(clSize, &val->invPerm));
3242: for (PetscInt i = 0; i < clSize; ++i) val->invPerm[clPerm[i]] = i;
3243: PetscFunctionReturn(PETSC_SUCCESS);
3244: }
3246: /*@
3247: PetscSectionSetClosurePermutation - Set the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3249: Not Collective
3251: Input Parameters:
3252: + section - The `PetscSection`
3253: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3254: . depth - Depth of points on which to apply the given permutation
3255: - perm - Permutation of the cell dof closure
3257: Level: intermediate
3259: Notes:
3260: The specified permutation will only be applied to points at depth whose closure size matches the length of perm. In a
3261: mixed-topology or variable-degree finite element space, this function can be called multiple times at each depth for
3262: each topology and degree.
3264: This approach assumes that (depth, len(perm)) uniquely identifies the desired permutation; this might not be true for
3265: exotic/enriched spaces on mixed topology meshes.
3267: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionGetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `DMPlexCreateClosureIndex()`, `PetscCopyMode`
3268: @*/
3269: PetscErrorCode PetscSectionSetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, IS perm)
3270: {
3271: const PetscInt *clPerm = NULL;
3272: PetscInt clSize = 0;
3274: PetscFunctionBegin;
3275: if (perm) {
3276: PetscCall(ISGetLocalSize(perm, &clSize));
3277: PetscCall(ISGetIndices(perm, &clPerm));
3278: }
3279: PetscCall(PetscSectionSetClosurePermutation_Internal(section, obj, depth, clSize, PETSC_COPY_VALUES, (PetscInt *)clPerm));
3280: if (perm) PetscCall(ISRestoreIndices(perm, &clPerm));
3281: PetscFunctionReturn(PETSC_SUCCESS);
3282: }
3284: static PetscErrorCode PetscSectionGetClosurePermutation_Private(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3285: {
3286: PetscFunctionBegin;
3287: if (section->clObj == obj) {
3288: PetscSectionClosurePermKey k = {depth, size};
3289: PetscSectionClosurePermVal v;
3291: PetscCall(PetscClPermGet(section->clHash, k, &v));
3292: if (perm) *perm = v.perm;
3293: } else {
3294: if (perm) *perm = NULL;
3295: }
3296: PetscFunctionReturn(PETSC_SUCCESS);
3297: }
3299: /*@
3300: PetscSectionGetClosurePermutation - Get the dof permutation for the closure of each cell in the section, meaning clPerm[newIndex] = oldIndex.
3302: Not Collective
3304: Input Parameters:
3305: + section - The `PetscSection`
3306: . obj - A `PetscObject` which serves as the key for this index (usually a DM)
3307: . depth - Depth stratum on which to obtain closure permutation
3308: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3310: Output Parameter:
3311: . perm - The dof closure permutation
3313: Level: intermediate
3315: Note:
3316: The user must destroy the `IS` that is returned.
3318: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureInversePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3319: @*/
3320: PetscErrorCode PetscSectionGetClosurePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3321: {
3322: const PetscInt *clPerm = NULL;
3324: PetscFunctionBegin;
3325: PetscCall(PetscSectionGetClosurePermutation_Private(section, obj, depth, clSize, &clPerm));
3326: PetscCheck(clPerm, PetscObjectComm(obj), PETSC_ERR_ARG_WRONG, "There is no closure permutation associated with this object for depth %" PetscInt_FMT " of size %" PetscInt_FMT, depth, clSize);
3327: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3328: PetscFunctionReturn(PETSC_SUCCESS);
3329: }
3331: PetscErrorCode PetscSectionGetClosureInversePermutation_Internal(PetscSection section, PetscObject obj, PetscInt depth, PetscInt size, const PetscInt *perm[])
3332: {
3333: PetscFunctionBegin;
3334: if (section->clObj == obj && section->clHash) {
3335: PetscSectionClosurePermKey k = {depth, size};
3336: PetscSectionClosurePermVal v;
3337: PetscCall(PetscClPermGet(section->clHash, k, &v));
3338: if (perm) *perm = v.invPerm;
3339: } else {
3340: if (perm) *perm = NULL;
3341: }
3342: PetscFunctionReturn(PETSC_SUCCESS);
3343: }
3345: /*@
3346: PetscSectionGetClosureInversePermutation - Get the inverse dof permutation for the closure of each cell in the section, meaning clPerm[oldIndex] = newIndex.
3348: Not Collective
3350: Input Parameters:
3351: + section - The `PetscSection`
3352: . obj - A `PetscObject` which serves as the key for this index (usually a `DM`)
3353: . depth - Depth stratum on which to obtain closure permutation
3354: - clSize - Closure size to be permuted (e.g., may vary with element topology and degree)
3356: Output Parameter:
3357: . perm - The dof closure permutation
3359: Level: intermediate
3361: Note:
3362: The user must destroy the `IS` that is returned.
3364: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetClosurePermutation()`, `PetscSectionGetClosureIndex()`, `PetscSectionSetClosureIndex()`, `DMPlexCreateClosureIndex()`
3365: @*/
3366: PetscErrorCode PetscSectionGetClosureInversePermutation(PetscSection section, PetscObject obj, PetscInt depth, PetscInt clSize, IS *perm)
3367: {
3368: const PetscInt *clPerm = NULL;
3370: PetscFunctionBegin;
3371: PetscCall(PetscSectionGetClosureInversePermutation_Internal(section, obj, depth, clSize, &clPerm));
3372: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPerm, PETSC_USE_POINTER, perm));
3373: PetscFunctionReturn(PETSC_SUCCESS);
3374: }
3376: /*@
3377: PetscSectionGetField - Get the `PetscSection` associated with a single field
3379: Input Parameters:
3380: + s - The `PetscSection`
3381: - field - The field number
3383: Output Parameter:
3384: . subs - The `PetscSection` for the given field, note the chart of `subs` is not set
3386: Level: intermediate
3388: Note:
3389: Does not increase the reference count of the selected sub-section. There is no matching `PetscSectionRestoreField()`
3391: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `IS`, `PetscSectionSetNumFields()`
3392: @*/
3393: PetscErrorCode PetscSectionGetField(PetscSection s, PetscInt field, PetscSection *subs)
3394: {
3395: PetscFunctionBegin;
3397: PetscAssertPointer(subs, 3);
3398: PetscSectionCheckValidField(field, s->numFields);
3399: *subs = s->field[field];
3400: PetscFunctionReturn(PETSC_SUCCESS);
3401: }
3403: PetscClassId PETSC_SECTION_SYM_CLASSID;
3404: PetscFunctionList PetscSectionSymList = NULL;
3406: /*@
3407: PetscSectionSymCreate - Creates an empty `PetscSectionSym` object.
3409: Collective
3411: Input Parameter:
3412: . comm - the MPI communicator
3414: Output Parameter:
3415: . sym - pointer to the new set of symmetries
3417: Level: developer
3419: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSectionSym`, `PetscSectionSymDestroy()`
3420: @*/
3421: PetscErrorCode PetscSectionSymCreate(MPI_Comm comm, PetscSectionSym *sym)
3422: {
3423: PetscFunctionBegin;
3424: PetscAssertPointer(sym, 2);
3425: PetscCall(ISInitializePackage());
3427: PetscCall(PetscHeaderCreate(*sym, PETSC_SECTION_SYM_CLASSID, "PetscSectionSym", "Section Symmetry", "IS", comm, PetscSectionSymDestroy, PetscSectionSymView));
3428: PetscFunctionReturn(PETSC_SUCCESS);
3429: }
3431: /*@
3432: PetscSectionSymSetType - Builds a `PetscSectionSym`, for a particular implementation.
3434: Collective
3436: Input Parameters:
3437: + sym - The section symmetry object
3438: - method - The name of the section symmetry type
3440: Level: developer
3442: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymGetType()`, `PetscSectionSymCreate()`
3443: @*/
3444: PetscErrorCode PetscSectionSymSetType(PetscSectionSym sym, PetscSectionSymType method)
3445: {
3446: PetscErrorCode (*r)(PetscSectionSym);
3447: PetscBool match;
3449: PetscFunctionBegin;
3451: PetscCall(PetscObjectTypeCompare((PetscObject)sym, method, &match));
3452: if (match) PetscFunctionReturn(PETSC_SUCCESS);
3454: PetscCall(PetscFunctionListFind(PetscSectionSymList, method, &r));
3455: PetscCheck(r, PetscObjectComm((PetscObject)sym), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSectionSym type: %s", method);
3456: PetscTryTypeMethod(sym, destroy);
3457: sym->ops->destroy = NULL;
3459: PetscCall((*r)(sym));
3460: PetscCall(PetscObjectChangeTypeName((PetscObject)sym, method));
3461: PetscFunctionReturn(PETSC_SUCCESS);
3462: }
3464: /*@
3465: PetscSectionSymGetType - Gets the section symmetry type name (as a string) from the `PetscSectionSym`.
3467: Not Collective
3469: Input Parameter:
3470: . sym - The section symmetry
3472: Output Parameter:
3473: . type - The index set type name
3475: Level: developer
3477: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymSetType()`, `PetscSectionSymCreate()`
3478: @*/
3479: PetscErrorCode PetscSectionSymGetType(PetscSectionSym sym, PetscSectionSymType *type)
3480: {
3481: PetscFunctionBegin;
3483: PetscAssertPointer(type, 2);
3484: *type = ((PetscObject)sym)->type_name;
3485: PetscFunctionReturn(PETSC_SUCCESS);
3486: }
3488: /*@C
3489: PetscSectionSymRegister - Registers a new section symmetry implementation
3491: Not Collective, No Fortran Support
3493: Input Parameters:
3494: + sname - The name of a new user-defined creation routine
3495: - function - The creation routine itself
3497: Level: developer
3499: Notes:
3500: `PetscSectionSymRegister()` may be called multiple times to add several user-defined vectors
3502: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymType`, `PetscSectionSymCreate()`, `PetscSectionSymSetType()`
3503: @*/
3504: PetscErrorCode PetscSectionSymRegister(const char sname[], PetscErrorCode (*function)(PetscSectionSym))
3505: {
3506: PetscFunctionBegin;
3507: PetscCall(ISInitializePackage());
3508: PetscCall(PetscFunctionListAdd(&PetscSectionSymList, sname, function));
3509: PetscFunctionReturn(PETSC_SUCCESS);
3510: }
3512: /*@
3513: PetscSectionSymDestroy - Destroys a section symmetry.
3515: Collective
3517: Input Parameter:
3518: . sym - the section symmetry
3520: Level: developer
3522: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`
3523: @*/
3524: PetscErrorCode PetscSectionSymDestroy(PetscSectionSym *sym)
3525: {
3526: SymWorkLink link, next;
3528: PetscFunctionBegin;
3529: if (!*sym) PetscFunctionReturn(PETSC_SUCCESS);
3531: if (--((PetscObject)*sym)->refct > 0) {
3532: *sym = NULL;
3533: PetscFunctionReturn(PETSC_SUCCESS);
3534: }
3535: PetscTryTypeMethod(*sym, destroy);
3536: PetscCheck(!(*sym)->workout, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Work array still checked out");
3537: for (link = (*sym)->workin; link; link = next) {
3538: PetscInt **perms = (PetscInt **)link->perms;
3539: PetscScalar **rots = (PetscScalar **)link->rots;
3540: PetscCall(PetscFree2(perms, rots));
3541: next = link->next;
3542: PetscCall(PetscFree(link));
3543: }
3544: (*sym)->workin = NULL;
3545: PetscCall(PetscHeaderDestroy(sym));
3546: PetscFunctionReturn(PETSC_SUCCESS);
3547: }
3549: /*@
3550: PetscSectionSymView - Displays a section symmetry
3552: Collective
3554: Input Parameters:
3555: + sym - the index set
3556: - viewer - viewer used to display the set, for example `PETSC_VIEWER_STDOUT_SELF`.
3558: Level: developer
3560: .seealso: `PetscSectionSym`, `PetscViewer`, `PetscViewerASCIIOpen()`
3561: @*/
3562: PetscErrorCode PetscSectionSymView(PetscSectionSym sym, PetscViewer viewer)
3563: {
3564: PetscFunctionBegin;
3566: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sym), &viewer));
3568: PetscCheckSameComm(sym, 1, viewer, 2);
3569: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sym, viewer));
3570: PetscTryTypeMethod(sym, view, viewer);
3571: PetscFunctionReturn(PETSC_SUCCESS);
3572: }
3574: /*@
3575: PetscSectionSetSym - Set the symmetries for the data referred to by the section
3577: Collective
3579: Input Parameters:
3580: + section - the section describing data layout
3581: - sym - the symmetry describing the affect of orientation on the access of the data
3583: Level: developer
3585: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetSym()`, `PetscSectionSymCreate()`
3586: @*/
3587: PetscErrorCode PetscSectionSetSym(PetscSection section, PetscSectionSym sym)
3588: {
3589: PetscFunctionBegin;
3591: PetscCall(PetscSectionSymDestroy(§ion->sym));
3592: if (sym) {
3594: PetscCheckSameComm(section, 1, sym, 2);
3595: PetscCall(PetscObjectReference((PetscObject)sym));
3596: }
3597: section->sym = sym;
3598: PetscFunctionReturn(PETSC_SUCCESS);
3599: }
3601: /*@
3602: PetscSectionGetSym - Get the symmetries for the data referred to by the section
3604: Not Collective
3606: Input Parameter:
3607: . section - the section describing data layout
3609: Output Parameter:
3610: . sym - the symmetry describing the affect of orientation on the access of the data, provided previously by `PetscSectionSetSym()`
3612: Level: developer
3614: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetSym()`, `PetscSectionSymCreate()`
3615: @*/
3616: PetscErrorCode PetscSectionGetSym(PetscSection section, PetscSectionSym *sym)
3617: {
3618: PetscFunctionBegin;
3620: *sym = section->sym;
3621: PetscFunctionReturn(PETSC_SUCCESS);
3622: }
3624: /*@
3625: PetscSectionSetFieldSym - Set the symmetries for the data referred to by a field of the section
3627: Collective
3629: Input Parameters:
3630: + section - the section describing data layout
3631: . field - the field number
3632: - sym - the symmetry describing the affect of orientation on the access of the data
3634: Level: developer
3636: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetFieldSym()`, `PetscSectionSymCreate()`
3637: @*/
3638: PetscErrorCode PetscSectionSetFieldSym(PetscSection section, PetscInt field, PetscSectionSym sym)
3639: {
3640: PetscFunctionBegin;
3642: PetscSectionCheckValidField(field, section->numFields);
3643: PetscCall(PetscSectionSetSym(section->field[field], sym));
3644: PetscFunctionReturn(PETSC_SUCCESS);
3645: }
3647: /*@
3648: PetscSectionGetFieldSym - Get the symmetries for the data referred to by a field of the section
3650: Collective
3652: Input Parameters:
3653: + section - the section describing data layout
3654: - field - the field number
3656: Output Parameter:
3657: . sym - the symmetry describing the affect of orientation on the access of the data
3659: Level: developer
3661: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetFieldSym()`, `PetscSectionSymCreate()`
3662: @*/
3663: PetscErrorCode PetscSectionGetFieldSym(PetscSection section, PetscInt field, PetscSectionSym *sym)
3664: {
3665: PetscFunctionBegin;
3667: PetscSectionCheckValidField(field, section->numFields);
3668: *sym = section->field[field]->sym;
3669: PetscFunctionReturn(PETSC_SUCCESS);
3670: }
3672: /*@C
3673: PetscSectionGetPointSyms - Get the symmetries for a set of points in a `PetscSection` under specific orientations.
3675: Not Collective
3677: Input Parameters:
3678: + section - the section
3679: . numPoints - the number of points
3680: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3681: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3682: context, see `DMPlexGetConeOrientation()`).
3684: Output Parameters:
3685: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3686: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3687: identity).
3689: Example of usage, gathering dofs into a local array (lArray) from a section array (sArray):
3690: .vb
3691: const PetscInt **perms;
3692: const PetscScalar **rots;
3693: PetscInt lOffset;
3695: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3696: for (i = 0, lOffset = 0; i < numPoints; i++) {
3697: PetscInt point = points[2*i], dof, sOffset;
3698: const PetscInt *perm = perms ? perms[i] : NULL;
3699: const PetscScalar *rot = rots ? rots[i] : NULL;
3701: PetscSectionGetDof(section,point,&dof);
3702: PetscSectionGetOffset(section,point,&sOffset);
3704: if (perm) { for (j = 0; j < dof; j++) lArray[lOffset + perm[j]] = sArray[sOffset + j]; }
3705: else { for (j = 0; j < dof; j++) lArray[lOffset + j ] = sArray[sOffset + j]; }
3706: if (rot) { for (j = 0; j < dof; j++) lArray[lOffset + j ] *= rot[j]; }
3707: lOffset += dof;
3708: }
3709: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3710: .ve
3712: Example of usage, adding dofs into a section array (sArray) from a local array (lArray):
3713: .vb
3714: const PetscInt **perms;
3715: const PetscScalar **rots;
3716: PetscInt lOffset;
3718: PetscSectionGetPointSyms(section,numPoints,points,&perms,&rots);
3719: for (i = 0, lOffset = 0; i < numPoints; i++) {
3720: PetscInt point = points[2*i], dof, sOffset;
3721: const PetscInt *perm = perms ? perms[i] : NULL;
3722: const PetscScalar *rot = rots ? rots[i] : NULL;
3724: PetscSectionGetDof(section,point,&dof);
3725: PetscSectionGetOffset(section,point,&sOff);
3727: if (perm) { for (j = 0; j < dof; j++) sArray[sOffset + j] += lArray[lOffset + perm[j]] * (rot ? PetscConj(rot[perm[j]]) : 1.); }
3728: else { for (j = 0; j < dof; j++) sArray[sOffset + j] += lArray[lOffset + j ] * (rot ? PetscConj(rot[ j ]) : 1.); }
3729: offset += dof;
3730: }
3731: PetscSectionRestorePointSyms(section,numPoints,points,&perms,&rots);
3732: .ve
3734: Level: developer
3736: Notes:
3737: `PetscSectionSetSym()` must have been previously called to provide the symmetries to the `PetscSection`
3739: Use `PetscSectionRestorePointSyms()` when finished with the data
3741: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3742: @*/
3743: PetscErrorCode PetscSectionGetPointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3744: {
3745: PetscSectionSym sym;
3747: PetscFunctionBegin;
3749: if (numPoints) PetscAssertPointer(points, 3);
3750: if (perms) *perms = NULL;
3751: if (rots) *rots = NULL;
3752: sym = section->sym;
3753: if (sym && (perms || rots)) {
3754: SymWorkLink link;
3756: if (sym->workin) {
3757: link = sym->workin;
3758: sym->workin = sym->workin->next;
3759: } else {
3760: PetscCall(PetscNew(&link));
3761: }
3762: if (numPoints > link->numPoints) {
3763: PetscInt **perms = (PetscInt **)link->perms;
3764: PetscScalar **rots = (PetscScalar **)link->rots;
3765: PetscCall(PetscFree2(perms, rots));
3766: PetscCall(PetscMalloc2(numPoints, (PetscInt ***)&link->perms, numPoints, (PetscScalar ***)&link->rots));
3767: link->numPoints = numPoints;
3768: }
3769: link->next = sym->workout;
3770: sym->workout = link;
3771: PetscCall(PetscArrayzero((PetscInt **)link->perms, numPoints));
3772: PetscCall(PetscArrayzero((PetscInt **)link->rots, numPoints));
3773: PetscUseTypeMethod(sym, getpoints, section, numPoints, points, link->perms, link->rots);
3774: if (perms) *perms = link->perms;
3775: if (rots) *rots = link->rots;
3776: }
3777: PetscFunctionReturn(PETSC_SUCCESS);
3778: }
3780: /*@C
3781: PetscSectionRestorePointSyms - Restore the symmetries returned by `PetscSectionGetPointSyms()`
3783: Not Collective
3785: Input Parameters:
3786: + section - the section
3787: . numPoints - the number of points
3788: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3789: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3790: context, see `DMPlexGetConeOrientation()`).
3791: . perms - The permutations for the given orientations: set to `NULL` at conclusion
3792: - rots - The field rotations symmetries for the given orientations: set to `NULL` at conclusion
3794: Level: developer
3796: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3797: @*/
3798: PetscErrorCode PetscSectionRestorePointSyms(PetscSection section, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3799: {
3800: PetscSectionSym sym;
3802: PetscFunctionBegin;
3804: sym = section->sym;
3805: if (sym && (perms || rots)) {
3806: SymWorkLink *p, link;
3808: for (p = &sym->workout; (link = *p); p = &link->next) {
3809: if ((perms && link->perms == *perms) || (rots && link->rots == *rots)) {
3810: *p = link->next;
3811: link->next = sym->workin;
3812: sym->workin = link;
3813: if (perms) *perms = NULL;
3814: if (rots) *rots = NULL;
3815: PetscFunctionReturn(PETSC_SUCCESS);
3816: }
3817: }
3818: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Array was not checked out");
3819: }
3820: PetscFunctionReturn(PETSC_SUCCESS);
3821: }
3823: /*@C
3824: PetscSectionGetFieldPointSyms - Get the symmetries for a set of points in a field of a `PetscSection` under specific orientations.
3826: Not Collective
3828: Input Parameters:
3829: + section - the section
3830: . field - the field of the section
3831: . numPoints - the number of points
3832: - points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3833: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3834: context, see `DMPlexGetConeOrientation()`).
3836: Output Parameters:
3837: + perms - The permutations for the given orientations (or `NULL` if there is no symmetry or the permutation is the identity).
3838: - rots - The field rotations symmetries for the given orientations (or `NULL` if there is no symmetry or the rotations are all
3839: identity).
3841: Level: developer
3843: Notes:
3844: `PetscSectionSetFieldSym()` must have been previously called to provide the symmetries to the `PetscSection`
3846: Use `PetscSectionRestoreFieldPointSyms()` when finished with the data
3848: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetPointSyms()`, `PetscSectionRestoreFieldPointSyms()`
3849: @*/
3850: PetscErrorCode PetscSectionGetFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3851: {
3852: PetscFunctionBegin;
3854: 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);
3855: PetscCall(PetscSectionGetPointSyms(section->field[field], numPoints, points, perms, rots));
3856: PetscFunctionReturn(PETSC_SUCCESS);
3857: }
3859: /*@C
3860: PetscSectionRestoreFieldPointSyms - Restore the symmetries returned by `PetscSectionGetFieldPointSyms()`
3862: Not Collective
3864: Input Parameters:
3865: + section - the section
3866: . field - the field number
3867: . numPoints - the number of points
3868: . points - an array of size 2 * `numPoints`, containing a list of (point, orientation) pairs. (An orientation is an
3869: arbitrary integer: its interpretation is up to sym. Orientations are used by `DM`: for their interpretation in that
3870: context, see `DMPlexGetConeOrientation()`).
3871: . perms - The permutations for the given orientations: set to NULL at conclusion
3872: - rots - The field rotations symmetries for the given orientations: set to NULL at conclusion
3874: Level: developer
3876: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionRestorePointSyms()`, `petscSectionGetFieldPointSyms()`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`
3877: @*/
3878: PetscErrorCode PetscSectionRestoreFieldPointSyms(PetscSection section, PetscInt field, PetscInt numPoints, const PetscInt *points, const PetscInt ***perms, const PetscScalar ***rots)
3879: {
3880: PetscFunctionBegin;
3882: 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);
3883: PetscCall(PetscSectionRestorePointSyms(section->field[field], numPoints, points, perms, rots));
3884: PetscFunctionReturn(PETSC_SUCCESS);
3885: }
3887: /*@
3888: PetscSectionSymCopy - Copy the symmetries, assuming that the point structure is compatible
3890: Not Collective
3892: Input Parameter:
3893: . sym - the `PetscSectionSym`
3895: Output Parameter:
3896: . nsym - the equivalent symmetries
3898: Level: developer
3900: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3901: @*/
3902: PetscErrorCode PetscSectionSymCopy(PetscSectionSym sym, PetscSectionSym nsym)
3903: {
3904: PetscFunctionBegin;
3907: PetscTryTypeMethod(sym, copy, nsym);
3908: PetscFunctionReturn(PETSC_SUCCESS);
3909: }
3911: /*@
3912: PetscSectionSymDistribute - Distribute the symmetries in accordance with the input `PetscSF`
3914: Collective
3916: Input Parameters:
3917: + sym - the `PetscSectionSym`
3918: - migrationSF - the distribution map from roots to leaves
3920: Output Parameter:
3921: . dsym - the redistributed symmetries
3923: Level: developer
3925: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSymCreate()`, `PetscSectionSetSym()`, `PetscSectionGetSym()`, `PetscSectionSymLabelSetStratum()`, `PetscSectionGetPointSyms()`
3926: @*/
3927: PetscErrorCode PetscSectionSymDistribute(PetscSectionSym sym, PetscSF migrationSF, PetscSectionSym *dsym)
3928: {
3929: PetscFunctionBegin;
3932: PetscAssertPointer(dsym, 3);
3933: PetscTryTypeMethod(sym, distribute, migrationSF, dsym);
3934: PetscFunctionReturn(PETSC_SUCCESS);
3935: }
3937: /*@
3938: PetscSectionGetUseFieldOffsets - Get the flag indicating if field offsets are used directly in a global section, rather than just the point offset
3940: Not Collective
3942: Input Parameter:
3943: . s - the global `PetscSection`
3945: Output Parameter:
3946: . flg - the flag
3948: Level: developer
3950: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3951: @*/
3952: PetscErrorCode PetscSectionGetUseFieldOffsets(PetscSection s, PetscBool *flg)
3953: {
3954: PetscFunctionBegin;
3956: *flg = s->useFieldOff;
3957: PetscFunctionReturn(PETSC_SUCCESS);
3958: }
3960: /*@
3961: PetscSectionSetUseFieldOffsets - Set the flag to use field offsets directly in a global section, rather than just the point offset
3963: Not Collective
3965: Input Parameters:
3966: + s - the global `PetscSection`
3967: - flg - the flag
3969: Level: developer
3971: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetUseFieldOffsets()`, `PetscSectionSetChart()`, `PetscSectionCreate()`
3972: @*/
3973: PetscErrorCode PetscSectionSetUseFieldOffsets(PetscSection s, PetscBool flg)
3974: {
3975: PetscFunctionBegin;
3977: s->useFieldOff = flg;
3978: PetscFunctionReturn(PETSC_SUCCESS);
3979: }
3981: #define PetscSectionExpandPoints_Loop(TYPE) \
3982: do { \
3983: PetscInt i, n, o0, o1, size; \
3984: TYPE *a0 = (TYPE *)origArray, *a1; \
3985: PetscCall(PetscSectionGetStorageSize(s, &size)); \
3986: PetscCall(PetscMalloc1(size, &a1)); \
3987: for (i = 0; i < npoints; i++) { \
3988: PetscCall(PetscSectionGetOffset(origSection, points_[i], &o0)); \
3989: PetscCall(PetscSectionGetOffset(s, i, &o1)); \
3990: PetscCall(PetscSectionGetDof(s, i, &n)); \
3991: PetscCall(PetscMemcpy(&a1[o1], &a0[o0], n * unitsize)); \
3992: } \
3993: *newArray = (void *)a1; \
3994: } while (0)
3996: /*@
3997: PetscSectionExtractDofsFromArray - Extracts elements of an array corresponding to DOFs of specified points.
3999: Not Collective
4001: Input Parameters:
4002: + origSection - the `PetscSection` describing the layout of the array
4003: . dataType - `MPI_Datatype` describing the data type of the array (currently only `MPIU_INT`, `MPIU_SCALAR`, `MPIU_REAL`)
4004: . origArray - the array; its size must be equal to the storage size of `origSection`
4005: - points - `IS` with points to extract; its indices must lie in the chart of `origSection`
4007: Output Parameters:
4008: + newSection - the new `PetscSection` describing the layout of the new array (with points renumbered 0,1,... but preserving numbers of DOFs)
4009: - newArray - the array of the extracted DOFs; its size is the storage size of `newSection`
4011: Level: developer
4013: .seealso: [PetscSection](ch_petscsection), `PetscSectionSym`, `PetscSectionGetChart()`, `PetscSectionGetDof()`, `PetscSectionGetStorageSize()`, `PetscSectionCreate()`
4014: @*/
4015: PetscErrorCode PetscSectionExtractDofsFromArray(PetscSection origSection, MPI_Datatype dataType, const void *origArray, IS points, PetscSection *newSection, void *newArray[])
4016: {
4017: PetscSection s;
4018: const PetscInt *points_;
4019: PetscInt i, n, npoints, pStart, pEnd;
4020: PetscMPIInt unitsize;
4022: PetscFunctionBegin;
4024: PetscAssertPointer(origArray, 3);
4026: if (newSection) PetscAssertPointer(newSection, 5);
4027: if (newArray) PetscAssertPointer(newArray, 6);
4028: PetscCallMPI(MPI_Type_size(dataType, &unitsize));
4029: PetscCall(ISGetLocalSize(points, &npoints));
4030: PetscCall(ISGetIndices(points, &points_));
4031: PetscCall(PetscSectionGetChart(origSection, &pStart, &pEnd));
4032: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &s));
4033: PetscCall(PetscSectionSetChart(s, 0, npoints));
4034: for (i = 0; i < npoints; i++) {
4035: 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);
4036: PetscCall(PetscSectionGetDof(origSection, points_[i], &n));
4037: PetscCall(PetscSectionSetDof(s, i, n));
4038: }
4039: PetscCall(PetscSectionSetUp(s));
4040: if (newArray) {
4041: if (dataType == MPIU_INT) {
4042: PetscSectionExpandPoints_Loop(PetscInt);
4043: } else if (dataType == MPIU_SCALAR) {
4044: PetscSectionExpandPoints_Loop(PetscScalar);
4045: } else if (dataType == MPIU_REAL) {
4046: PetscSectionExpandPoints_Loop(PetscReal);
4047: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for this MPI_Datatype");
4048: }
4049: if (newSection) {
4050: *newSection = s;
4051: } else {
4052: PetscCall(PetscSectionDestroy(&s));
4053: }
4054: PetscCall(ISRestoreIndices(points, &points_));
4055: PetscFunctionReturn(PETSC_SUCCESS);
4056: }
4058: /*@
4059: PetscSectionMigrateData - Migrate data described by a `PetscSection` using a `PetscSF` that defines a original-to-new (root-to-leaf) point mapping
4061: Collective
4063: Input Parameters:
4064: + migratePointSF - defines the mapping (communication) of the root points to the leaf points
4065: . datatype - the type of data
4066: . rootSection - the `PetscSection` that describes the data layout on the root points (how many dof and what fields are associated with each root point)
4067: - rootData - the existing data array described by `rootSection`, may be `NULL` is storage size of `rootSection` is zero
4069: Output Parameters:
4070: + leafSection - the new `PetscSection` that describes the data layout on the leaf points
4071: . leafData - the redistributed data array that is associated with the leaf points
4072: - migrateDataSF - defines the mapping (communication) of the `rootData` array to the `leafData` array, may be `NULL` if not needed
4074: Level: advanced
4076: Notes:
4077: This function can best be thought of as applying `PetscSFBcastBegin()` to an array described by a `PetscSection`.
4078: While `PetscSFBcastBegin()` is limited to broadcasting data that is of the same size for every index, this function allows the data to be a different size for each index.
4079: The size and layout of that irregularly sized data before and after `PetscSFBcastBegin()` is described by the `rootSection` and `leafSection`, respectively.
4081: This function combines `PetscSFDistributeSection()`, `PetscSFCreateSectionSF()`, and `PetscSFBcastBegin()`/`PetscSFBcastEnd()` into a single call.
4082: `migrateDataSF` can be used to repeat the `PetscSFBcastBegin()`/`PetscSFBcastEnd()` on a different data array described by the same `rootSection`.
4084: This should not be used for global-to-local type communication patterns.
4085: For this use case, see `PetscSectionCreateGlobalSection()` and `PetscSFSetGraphSection()`.
4087: .seealso: [PetscSection](ch_petscsection), `PetscSection`, `PetscSFDistributeSection()`, `PetscSFCreateSectionSF()`, `DMPlexDistributeData()`
4088: @*/
4089: PetscErrorCode PetscSectionMigrateData(PetscSF migratePointSF, MPI_Datatype datatype, PetscSection rootSection, const void *rootData, PetscSection leafSection, void *leafData[], PetscSF *migrateDataSF)
4090: {
4091: PetscSF fieldSF;
4092: PetscInt *remoteOffsets, fieldSize;
4093: PetscMPIInt dataSize;
4095: PetscFunctionBegin;
4098: if (rootData) PetscAssertPointer(rootData, 4);
4099: else {
4100: PetscInt size;
4101: PetscCall(PetscSectionGetStorageSize(rootSection, &size));
4102: PetscCheck(size == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "originalData may be NULL iff the storage size of originalSection is zero, but is %" PetscInt_FMT, size);
4103: }
4105: PetscAssertPointer(leafData, 6);
4106: if (migrateDataSF) PetscAssertPointer(migrateDataSF, 7);
4108: PetscCall(PetscSFDistributeSection(migratePointSF, rootSection, &remoteOffsets, leafSection));
4109: PetscCall(PetscSFCreateSectionSF(migratePointSF, rootSection, remoteOffsets, leafSection, &fieldSF));
4110: PetscCall(PetscFree(remoteOffsets));
4112: PetscCall(PetscSectionGetStorageSize(leafSection, &fieldSize));
4113: PetscCallMPI(MPI_Type_size(datatype, &dataSize));
4114: PetscCall(PetscMalloc(fieldSize * dataSize, leafData));
4115: PetscCall(PetscSFBcastBegin(fieldSF, datatype, rootData, *leafData, MPI_REPLACE));
4116: PetscCall(PetscSFBcastEnd(fieldSF, datatype, rootData, *leafData, MPI_REPLACE));
4118: if (migrateDataSF) *migrateDataSF = fieldSF;
4119: else PetscCall(PetscSFDestroy(&fieldSF));
4120: PetscFunctionReturn(PETSC_SUCCESS);
4121: }